function get_euvi_overlap,hdr,soho=soho,negative=negative,ushelio=ushelio
;+
; $Id: get_euvi_overlap.pro,v 1.4 2011/08/09 18:46:40 mcnutt Exp $
;
; Project : STEREO SECCHI
;
; Name    : get_euvi_overlap
;
; Purpose : returns mask of overlap regions where overlap where(overlap eq 255) 
;
; Use     : IDL>overlap=get_euvi_overlap(hdr)
;
; Inputs  :
;    hdr = secchi image hdr
;
; Outputs :
;    masks array of (hdr.naxis1,hdr.naxis2)
;
; Keywords:
;    /soho : retrun overlap with so instead of other sc
;    /negative : returns a mask where overlap region is eq to 0 not 255
;    /ushelio  selects overlap using the carrington longitudes return by euvi_heliographic
;
; Common  : None
;
; Restrictions: only works for full disk images. may(?) work on partial images if you use the keyword ushelio
;
; Side effects:
;
; Category    : Image analysis
;
; Prev. Hist. :
;
; Written     : Lynn Simpson 
;               based on Karls dons_euvi and longitude calculations from ssw/soho/lasco/idl/display/sungrid.pro
;
; $Log: get_euvi_overlap.pro,v $
; Revision 1.4  2011/08/09 18:46:40  mcnutt
; corrected for days with no overlap
;
; Revision 1.3  2010/06/10 16:48:15  mcnutt
; includes file type in filename for sccfindfits so not to return cal images
;
; Revision 1.2  2009/10/06 18:53:50  mcnutt
; added ushelio option
;
; Revision 1.1  2009/06/18 14:23:42  mcnutt
; return mask of overlap region
;
;

;-----------------
; DONS_EUVI (will get a new name once released!)
; originaly Written by: Karl Battams, NRL/Interferometrics,Inc, Oct26,2007
;-----------------

; Give this routine the name of a SECCHI EUVI image and it will 
; display that image and indicate the sector of the solar disk 
; that is only visible to itself (i.e. the bit that's "round the
; corner" on the other spacecraft. 

; I need to modify it so that it somehow highlights the sector in
; question, but for the most part this is pretty slick. (I just 
; hope it's correct!!  EDIT: It's not...)

if keyword_set(soho) and keyword_set(ushelio) then begin
   ushelio=0
   print,'ushelio keyword does not workwith SOHO Yet!'
endif

IF ~keyword_set(ushelio) THEN BEGIN

wave=hdr.WAVELNTH


obs=strmid(hdr.OBSRVTRY,7,1)
if obs eq 'A' then other ='B' else other='A'
if keyword_set(soho) then begin
  other='SOHO' 
  load_soho_spice
endif
sepang=GET_STEREO_SEP_ANGLE( hdr.DATE_OBS, obs,other,/degrees )
print, "S/C Sep. Angle: "+strcompress(sepang)

if hdr.crln_obs lt 0 then hdr.crln_obs = hdr.crln_obs +360 ; correct headers through Feb 09 which have negative longitudes
IF hdr.OBSRVTRY EQ 'STEREO_B' THEN $
    longs=[hdr.crln_obs-90+sepang,hdr.crln_obs+90] ELSE $
    longs=[hdr.crln_obs-90,hdr.crln_obs+90-sepang]


ckl=where(longs gt 360) 
if ckl(0) gt -1 then longs(ckl)=longs(ckl)-360

ckl=where(longs lt 0) 
if ckl(0) gt -1 then longs(ckl)=longs(ckl)+360


;###############################################################################
; Define some constants for the grid
;###############################################################################

n=200               ; The normal number of points out of which a grid line consists
r=(1.02*( ((2*hdr.RSUN)/hdr.CDELT1) ))/2.  ;solar radious

; create the array with the latitudes and longitudes to be drawn 

B_arr=!DTOR*[-90, -75, -60, -45, -30, -15, 0, 15, 30, 45, 60, 75, 90]

L_arr=!DTOR*longs

n_b=N_ElEMENTS(B_arr)
n_L=N_ELEMENTS(L_arr)

;###############################################################################
; Erzeugung der Drehmatrix
;###############################################################################
crot=!DTOR*(hdr.CROTA)
Bz=!DTOR*hdr.CRLT_OBS
Lz=!DTOR*hdr.crln_obs
;stop

; Transformation to rotate the sun to a certain central meridian
D1=[[ cos(Lz), 0,-sin(Lz)], $
    [       0, 1,       0], $
    [sin(Lz), 0, cos(Lz)]]

; Transformation to take into account the central latitude
D2=[[ 1,      0,        0], $
    [0, cos(Bz), -sin(Bz)], $
    [0, sin(Bz),  cos(Bz)]]

; Transformation to take into account the position angle
D3=[[cos(crot), -sin(crot), 0], $
    [sin(crot),  cos(crot), 0], $
    [     0,       0, 1]]

; Final rotation of the coordinate points
D=D1#D2#D3



;###############################################################################
; Longitude circles
;###############################################################################

b=!PI*(FINDGEN(n+1))/n
b=b-!pi/2


olddev=!D.NAME
SET_PLOT,'z'
DEVICE,SET_RESOLUTION=[hdr.naxis1,hdr.naxis2]
ERASE

thick=1

 FOR i=0,n_l-1 DO BEGIN
    l=L_arr(i)
    IF ((!RADEG*l MOD 90) EQ 0) THEN lthick=2*thick ELSE lthick=thick
    
; determine the [x,y,z] coordinates for the point on the circle
    x=r*sin(l)*cos(b)
    y=r*sin(b)
    z=r*cos(l)*cos(b)
    
; rotate the points to reflect the (P,B0,L0) orientation (matrix
; multiplication with the rotation matrix)
    xn=D(0,0)*x+D(1,0)*y+D(2,0)*z
    yn=D(0,1)*x+D(1,1)*y+D(2,1)*z
    zn=D(0,2)*x+D(1,2)*y+D(2,2)*z

    xn=xn+hdr.CRPIX1
    yn=yn+hdr.CRPIX2
    
    PLOTS,xn,yn,/DEVICE,THICK=lthick, COLOR=255
 ENDFOR

gd=TVRD()
SET_PLOT,olddev

;tv,gd
nonzero=where(gd EQ 255,count)
arr_ind=array_indices(gd,nonzero)

; This bit trims off loose ends of lines via the principal that a no-zero pixel is 
; part of a loose end if it does not neighbor at least two adjacent non-zero pixels
; I've commented it out because it relies on an entirely enclosed area, otherwise it 
; just erases whatever boundary you've defined!  It's not necessary for this app
; anyway, but once I wrote it, I decided to keep it hanging around. You never know...
    
;changecount=1
;while (changecount GT 0) DO BEGIN
;changecount=0
;for xv=1,outsize-2 DO BEGIN
;    for yv=1,outsize-2 DO BEGIN
;    	pix=gd[xv,yv]
;	IF pix GT 0 THEN BEGIN
;	    tot=total(gd[xv-1:xv+1,yv-1:yv+1])
;	    IF tot LE (2*255) THEN BEGIN
;	    	gd[xv,yv]=0
;		changecount+=1
;	    ENDIF
;	ENDIF
;    ENDFOR
;ENDFOR
;print,changecount
;ENDWHILE

; Here, we are going to fill in the boundary (ROI) with '255'.  I can't find an existing IDL routine
; that does exactly what I want, so the following works nicely.
; First it crops the array down to a box that tightly bounds our ROI.  Then it considers each 
; individual pixel within the crop and determines if it's inside the ROI.  How does it do this?
; Well, it assumes that for any point inside a ROI, if you traverse the row or column occupied by
; that pixel in all four directions, you will hit a bound at some point.  If that's true, the point 
; must be inside the ROI.  This works as long as the ROI is continuous and there are no other features
; in the image

xmax=max(arr_ind[0,*]+1)
xmin=min(arr_ind[0,*]-1)
ymax=max(arr_ind[1,*]+1)
ymin=min(arr_ind[1,*]-1)
crop=gd[xmin:xmax,ymin:ymax]   ; get the crop

sz=size(crop)
copy=crop   ; make a copy of the crop
cx=sz[1]
cy=sz[2]

FOR xv=0,cx-1 DO BEGIN
  ck=where(crop(xv,*) gt 0)
  if ck(0) gt -1 then copy(xv,min(ck):max(ck))=255
ENDFOR
FOR yv=0,cy-1 DO BEGIN
  ck=where(crop(*,yv) gt 0)
  if ck(0) gt -1 then copy(min(ck):max(ck),yv)=255
ENDFOR

;FOR xv=0,cx-1 DO BEGIN
;    FOR yv=0,cy-1 DO BEGIN
;    	pix=crop[xv,yv]     	    ; consider the pixel
;	IF (pix NE 255) THEN BEGIN
;	    uprow=crop[xv:*,yv]     ; get the rows/cols in all directions from the pixel
;	    downrow=crop[0:xv,yv]
;	    upcol=crop[xv,0:yv]
;	    downcol=crop[xv,yv:*]
;	    colsum=float(max(uprow))+float(max(downrow))+float(max(downcol))+float(max(upcol))
	    ; if the sum of those columns GE 1020 (255*4) then it's an internal point
;	    IF (colsum GE 1020) THEN copy[xv,yv]=255 
;	ENDIF
;    ENDFOR
;ENDFOR
;stop

gd[xmin:xmax,ymin:ymax]=copy
; Take out the next four lines to 'reverse' the mask overlay
ENDIF ELSE BEGIN
  file=hdr.filename
  if strpos(file,'A.fts') gt -1 then osc ='B' else osc='A'
  strput,file,osc,strpos(file,'.f')-1
  ;strput,file,'?',strpos(file,'.f')-5
  strput,file,'??????',strpos(file,'.f')-12
  fname = sccfindfits(file)

  if fname ne '' then begin
   jk=sccreadfits(fname,scch,/nodata)
   wcs2=fitshead2wcs(scch)
  ENDIF else wcs2= future_wcshdr(osc,hdr.detector,strmid(hdr.date_obs,0,10),strmid(hdr.date_obs,11,strlen(hdr.date_obs)))

  wcs1=fitshead2wcs(hdr)
  
  clt1=euvi_heliographic(wcs1,/CARRINGTON)                                  
  pxs=intarr(2,wcs2.naxis(0))
  pixs=lonarr(2,wcs2.naxis(0))
  pixs(1,*)=intarr(wcs2.naxis(1))+round(wcs2.crpix(1))
  pixs(0,*)=indgen(wcs2.naxis(0))
  clt2=euvi_heliographic(wcs2,pixs,/CARRINGTON)                                  

  cllng1=reform(clt1(0,*,*))
  cllng1(where (~Finite(cllng1)))=-3000

  cllng2=reform(clt2(0,*))
  cllng2(where (~Finite(cllng2)))=-3000
  side=cllng2(where(cllng2 gt -1))
  limb=[side(0),side(n_Elements(side)-1)]

  gd=bytarr(wcs1.naxis(0),wcs1.naxis(1))
  if limb(0) gt limb(1) then begin
     cllng1(where(cllng1 lt 180))=cllng1(where(cllng1 lt 180))+360
     limb(1)=limb(1)+360
  endif

  tst=where(cllng1 ge limb(0) and cllng1 le limb(1),cnt)
  if cnt ge 1 then gd(tst)=255

ENDELSE


if keyword_set(negative) then begin
  white=where(gd EQ 255,cntw)
  black = where(gd LT 255,cntb)
  if cntw ge 1 then gd[white]=0
  if cntb ge 1 then gd[black]=255
endif

return,gd

END