;+
function rt_carrmapmaker, yymmdd, num_r, rad, limb, wlimb, hdr, disp, saveset, nextrot

;
;PURPOSE: 
;	Restores partially finished Carrington Maps for a given limb/camera 
;	and appends the passed day's worth of images.
;
;INPUTS:
;	yymmdd	STRING:	date
;	num_r	INT:	number of radii to make maps for-passed empty
;	rad	FLTARR(num_r):		radii for maps-passed empty
;	cmap	FLTARR(mapsize,181,num_r):	passed undefined
;	limb	INTARR:	0 or 1 tells which limb for carrdate procedure
;	disp	INT:	1 to display images
;	wlimb	STRING:	'wl' or 'el'
;	hdr	STRUCT: dummy header sent back to carrmap3 for fits header
;	saveset	STRING:	filename of saved set of carrmapmaker
;	nextrot	INT	1 if going on to the next rotation
;
;
; OUTPUTS:	rad = fltarr(num_r)	:	radii of maps
;		num_r	: number of radii
;		RESULT = fltarr(mapsize,181,num_r):	carrington maps, one per radius
;		
;
;ROUTINES CALLED:
;    carrdate2.pro
;    getc2c3norm2.pro
;    savestrips.pro
;
; AUTHOR:	Nathan Rich, NRL, Nov. 1996
;		Julia Kraemer, NRL, June 7, 1996
;
; MODIFIED:
;	11/12/98 NBR	modified from carrmapmaker2.pro
;	12/09/98 NBR	do cn=cn+1 in beginning
;	04/16/99 NBR	make into a function
;	06/08/99 NBR	Put check for m LE 0 (median)
;	11/23/99 NBR	Change strip end skip criteria
;	02/2000  NBR	GET_BKG=2 for any_year
;	12/17/01 NBR	Use GET_BKG=1 for current year
;
; 06/08/07, @(#)rt_carrmapmaker.pro	1.3
;
;-

COMMON warpit, carrmap_x, stripn0, time0, duration, mapsize,  rat, startime, endtime, normalize, sta_radius, end_radius, sta_center, end_center, rollcor, platescl, cn
COMMON map_display, mapmax,mapmin

rotangle = 0
df=2
disp_factor=df
stripends = fltarr(22)
end_center=fltarr(2)

restore, saveset
IF nextrot THEN cn = cn+1
nextrot=0
qldir=getenv('QL_IMG')+'/level_05/'
filelist= findfile(qldir+yymmdd+'/'+STRLOWCASE(cam)+'/*.fts')
numim=n_elements(filelist)
if numim le 1 then BEGIN
	print,'One or zero images for ',cam,' on ',yymmdd,': Returning'
	nextrot = -1
	RETURN,-1
ENDIF
IF disp EQ 1 then BEGIN
	window,2,xsize=512,ysize=512
	window,0,xsize=mapsize,ysize=181
ENDIF

num_r=n_elements(rad)
time0s = carrdate2(cn,wl=limb(0), el=limb(1), cmp = limb(2), clong = limb(3))
time1s = carrdate2(cn+1,wl=limb(0), el=limb(1), cmp = limb(2), clong = limb(3))
startime=utc2str(time0s,/ecs)
endtime =utc2str(time1s,/ecs)
t_start = time0s.mjd + double(time0s.time/86400000.)
t_end = time1s.mjd + time1s.time/86400000.	
OPENW,l1,cam+wlimb+'_'+TRIM(STRING(cn))+'_'+TRIM(STRING(rad(0)))+'.log' ,/get_lun, /append
get_utc,now,/ecs
printf,l1, '*** rt_carrmapmaker resumed at ',now

IF camnum EQ 1 THEN deltend = 0.1
IF camnum EQ 2 THEN deltend = 0.008

FOR k = 0,numim-1 DO BEGIN
   print,'carrmapmaker: k = ',k
   toskip = 'n'
   filename = filelist(k)
   filen = filename
   ;IF (mapdelta EQ 0 AND carrmap_x NE 0) THEN BEGIN
;	print,'Skipping file ',filen
;	GOTO, skipped
   ;ENDIF
   print, 'Reading file ',filen  
   nothing = lasco_readfits(filen,hdr,/no_img)
   cam = hdr.detector
   IF (( cam EQ 'C3' and hdr.filter NE 'Clear') or $
       ( cam EQ 'C2' and hdr.filter NE 'Orange')) or $
      hdr.polar NE 'Clear' THEN GOTO, skipped
   obs_timen = hdr.date_obs+' '+hdr.time_obs

;
;  ** Begin carrwarp1. **
;
   dayn = str2utc(obs_timen)
   tn = DOUBLE(dayn.mjd + dayn.time/86400000.)
 				; time is milliseconds, changed to day units
   				; day is modified julian day
   tn0 = t_start + carrmap_x/pixperday
   deltaday = tn - tn0		; time between last column and present image
   mapdelta = ROUND(deltaday*pixperday)
   print, 'tn =',tn,'  deltaday =',deltaday,'  mapdelta =',mapdelta
;stop
   IF mapdelta LE 0 THEN BEGIN
	printf,l1,'Skipped ',filen,': Out of Order or too close to previous.'
	print,'Skipping ',filen,': Out of Order or too close to previous.'
	GOTO, skipped
   ENDIF

   IF tn GE t_end THEN BEGIN
	printf,l1,'Skipped ',filen,': In next rotation.'
	print,'Skipping ',filen,': In next rotation.'
	nextrot = 1
	GOTO, skipped
   ENDIF


;
;  ** End carrwarp1.   **
;
   im0 = getc123(filen,hdr,GET_BKG=1, RUNNING_DIFF=running_diff)
   ;im = reduce_std_size(im0,hdr,/full)		; sizes im to 1024x1024 
   im = im0

  ; number = strmid(filename,0,8)
  ; stripfile = findfile('/net/corona/export/raid2/ew_strips/' +yymmdd +'/' +cam +'/' +number +'_' +wlimb +'.fts')

   utcday = dayn
   yymmdd = utc2yymmdd(utcday)
   utcdate = dayn
   solar_ephem,yymmdd,radius=radius,/soho
   rconv = radius*3600/platescl
   ;IF hdr.naxis1 LT 700 THEN sclfactor=2 ELSE sclfactor=1
   cx = hdr.crpix1*hdr.lebxsum
   cy = hdr.crpix2*hdr.lebysum

   szim = SIZE(im)
   IF szim(0) EQ 0 THEN BEGIN
	print,'File not found--skipping'
	printf,l1,'File not found--skipping'
	GOTO, skipped
   ENDIF

   pnt = GET_SC_POINT(obs_timen)
   rotangle = -1*pnt.sc_roll	; degrees
   IF ABS(rotangle) GT 1 THEN BEGIN
	im =rot(temporary(im),rotangle,1,cx,cy,/pivot,/interp)
	IF carrmap_x LE 1 THEN rollcor(0)=rotangle
	print,'Rotating ',filename,rotangle,' degrees CW.'
	printf,l1,'Rotating ',filename,rotangle,' degrees CW.'
   ENDIF

   IF szim(1) GT 600 THEN m = median(im(*,y21:y12)) $
	ELSE m=median(im(*,y21/2:y12/2))
   mx = max(im0)
   help,m

   IF disp EQ 1 THEN BEGIN
	IF szim(1) GT 600 THEN BEGIN
	   xs = szim(1)/2
	   ys = szim(2)/2
	   disp_im = rebin(im,xs,ys)
	ENDIF ELSE disp_im = im
	disp_im(x11/df:x12/df, y11/df) = mx
	disp_im(x11/df:x12/df, y12/df) = mx
	disp_im(x21/df:x22/df,y21/df) = mx
	disp_im(x21/df:x22/df,y22/df) = mx
	wset,2
 ;  	IF pol EQ 'p' OR pol EQ 'u' THEN scale = 80
   	tvscl,disp_im>mapmin<mapmax
   ENDIF
   IF m LE 0 THEN BEGIN
	print,'Skipping ',filen,': median LE zero.'
	GOTO, skipped
   ENDIF

;
;   IF cam NE 'c1' THEN BEGIN
;   	IF (m LE 1.1*median0 AND m GE 0.9*median0) THEN median0 = m ELSE $
;   	IF (m LE 1.1*median1 AND m GE 0.9*median1) THEN median1 = m $
;	ELSE BEGIN
;	   skip='y'
;;	   read,'Skip this one? [y/n]',skip
;	   IF skip EQ 'y' THEN BEGIN
;	   	print,'Skipping file: Median out of range.'
;	   	printf,l1,'Skipping file : Median out of range.'
;		skiplast = skiplast + 1
;		wait,2
;	   	GOTO,skipped
;	   ENDIF 
;   	ENDELSE
;   ENDIF


;  ** Begin intensity normalization. **
   IF normalize EQ 'T' AND carrmap_x EQ 0 THEN BEGIN
	box_img = DOUBLE(im(x11:x12,y11:y12))
        ;good = WHERE(box_img NE 0 AND box_img LT 1.4)
        ;IF (good(0) GE 0) THEN box_ref=TOTAL(box_img(good))/N_ELEMENTS(good) ELSE box=0
	box_ref = median(box_img)
   ENDIF

   a1=1
   a2=1
int_corr=4.3
;im=(4.3/int_corr)*im 
   IF normalize EQ 'T' THEN BEGIN	
;   IF cam EQ 'nothing' THEN BEGIN
       squ1=DOUBLE(im(x11:x12,y11:y12))
       squ2=DOUBLE(im(x21:x22,y21:y22))
   
   	ch1=where(squ1 NE 0 AND squ1 LT 1.4,null1)
   	ch2=where(squ2 NE 0 AND squ2 LT 1.4,null2)
	IF (ch1(0) LT 0 OR ch2(0) LT 0) THEN BEGIN
		print,'Skipping file: box median out of range.'
		printf,l1,'Skipping file: box median out of range.'
		wait,2
		GOTO,skipped
	ENDIF
	IF null1 LT 2 THEN a1=0
	IF null2 LT 2 THEN a2=0
   	IF a1 EQ a2 THEN BEGIN
		a1 = 0.5
		a2 = 0.5
   	ENDIF
   	;avg1=TOTAL(squ1(ch1))/null1
   	;avg2=TOTAL(squ2(ch2))/null2
	avg1 = median(squ1)
	avg2 = median(squ2)
   	int_corr=a1*avg1 + a2*avg2
   	print,'int_corr = ',int_corr,'    toplft =',null1,'    botrt =',null2
	factor = box_ref/int_corr	;** box_ref 0.17 for c1?
   	IF int_corr GT 0 THEN im=factor*im ELSE BEGIN
		print,'Skipping file: box median out of range.'
		printf,l1,'Skipping file: box median out of range.'
		wait,2
		GOTO,skipped
   	ENDELSE		; factor is about im = 0.12*im for c2, 0.15*im for c3?
   ENDIF
	
;  ** End intensity normalization. **

   westlimb = strpos(wlimb,'wl')
   FOR r=0,num_r-1 DO BEGIN               ; cycles through each radius
	r0=rad(r)*rconv
	rn=r0
	
; 	** Take a strip for each radius. **

	FOR i=0,180 DO BEGIN
	   theta=i-90
	   th=theta*!pi/180.
	   IF westlimb NE -1 THEN xcoor=r0*cos(th)+cx ELSE $
	   	xcoor=-r0*cos(th)+cx	
		;** NBR, 4/15/98, east limb maps unless wlimb has 'wl' in it
	   xwcoor=cx+r0*cos(th)
	   xecoor=cx-r0*cos(th)
	   ycoor=r0*sin(th)+cy
           strp(r,i)=INTERPOLATE(im,xcoor,ycoor)
	   IF disp EQ 1 THEN BEGIN
	   	disp_im(xecoor/disp_factor,ycoor/disp_factor)=mx
	   	disp_im(xwcoor/disp_factor,ycoor/disp_factor)=mx
	   ENDIF
	ENDFOR   ;	i loop
	
; 	** End take strip. **
;
;	w = where(stripn0)
;	IF w(0) LT 0 THEN stripn0(r,*) = strp(r,*)
;	;maxmin,strip(r,*)
;	strp(r,*) = strp(r,*)*avg_strip_ends(r)/newmedend
;	strp(r,*) = avg_med(r)*(strp(r,*)/median(strp(r,*)))
;
   ENDFOR	; r loop

;   savestrips,strp,hdr,factor,rad,wlimb

;IF k LT 2 THEN stop
printf,l1,filename+':  carrmap_x ='+STRING(carrmap_x)+',  rconv='+STRING(rconv) +'  median(im) ='+STRING(m)+'  max(im) =' +STRING(mx) +'  int_corr =' +STRING(int_corr)
print,filename+':  carrmap_x ='+STRING(carrmap_x)+',  rconv='+STRING(rconv) +'  median(im) ='+STRING(m)+'  max(im) =' +STRING(mx)		

   IF disp EQ 1 THEN BEGIN
	wset,2
 ;  	IF pol EQ 'p' OR pol EQ 'u' THEN scale = 80
   	tvscl,disp_im>mapmin<mapmax
   ENDIF
   ;IF k LT 54 and k GT 44 THEN stop

print
;; ** BEgin  check strip **

   stripends(0:10)=strp(3,0:10)
   stripends(11:21)=strp(3,170:180)
   newmedend = median(stripends)
   newmed = median(strp(3,*))
   oldmed = median(stripn0(3,*))
   print,'ends =',newmedend,'   strip median =',newmed
   printf,l1,'strip median =',newmed
   	IF newmedend LT 0 OR newmed LT 0 THEN BEGIN
		print,'Strip or end is negative; skipping image.'
		wait,2
		printf,l1,'Strip or end is negative; skipping image.'
		skiplast = skiplast + 1
		GOTO, skipped
	ENDIF
w=where(oldmedends)
   ww =where(oldstripmeds)
   IF w(0) LT 0 THEN est_val = avg_strip_ends(camnum) ELSE BEGIN
   	y = oldmedends(w(0):19)
	sy = size(y)
   	wsz = size(w)
   	IF wsz(1) EQ 1 THEN est_val = y(0) ELSE BEGIN
   	   ;x = lindgen(wsz(1))
   	   ;fun = poly_fit(x,y,1,yfit,yband,sigma)
	   ;est_val = fun(0) + fun(1)*wsz(1)
	   est_val = total(y)/sy(1)
   	ENDELSE
   ENDELSE
   IF ww(0) LT 0 THEN str_est_val = newmed ELSE BEGIN
   	yy = oldstripmeds(ww(0):9)
	syy= size(yy)
   	wwsz = size(ww)
   	IF wwsz(1) EQ 1 THEN str_est_val = yy(0) ELSE BEGIN
   	   ;xx = lindgen(wwsz(1))
   	   ;fun2 = poly_fit(xx,yy,1,yyfit,yyband,ssigma)
	   ;str_est_val = fun2(0) + fun2(1)*wwsz(1)
	   str_est_val = total(yy)/syy(1)
   	ENDELSE
   ENDELSE
   print,'est_val =',est_val,',   str_est_val =',str_est_val
   wait,2
   ;IF carrmap_x EQ 160 THEN stop
   IF ( (newmed LT str_est_val*(1-deltend)) OR $
	(newmed GT str_est_val*(1+(deltend))) ) THEN BEGIN
;		read,' Skip this one? (Default is yes) ',toskip
;		IF toskip EQ '' THEN toskip = 'y'
	toskip='y'
   ENDIF
			; changed from 3*str_est_val, 6/8/98

   IF ( (newmedend LT est_val*(1-(deltend+0.05))) OR $
	(newmedend GT est_val*(1+(deltend+0.05))) ) THEN toskip = 'y'
			; changed from 3*est_val, 6/8/98
			; changed from +0.004 to +0.05, 11/16/99
   IF toskip EQ 'y' THEN Begin
	IF (carrmap_x LT 5) THEN toskip='n'
   	IF mapdelta GT 5 and ( (newmedend GT est_val*(1-(deltend+0.5))) and $
	(newmedend LT est_val*(1+(deltend+0.5))) ) THEN BEGIN
		; nbr, 11/23/99
		toskip = 'n'
		print,'Mapdelta GT 5; not skipping'
	ENDIF
	wait,2
   
	IF toskip NE 'n' THEN BEGIN
		skiplast = skiplast+1
		print,'Strip is too different; skipping image.'
		wait,2
		printf,l1,'Strip is too different; skipping image.'
		GOTO, skipped
	ENDIF
   ENDIF
   skiplast = 0

;  ** END check strip **

   carrmap_x = carrmap_x + mapdelta
   IF carrmap_x GT mapsize-1 THEN carrmap_x = mapsize -1

   ;IF k EQ 0  THEN stripn0 = strp
   FOR r=0,num_r-1 DO BEGIN               ; cycles through each radius
;	** Begin carrwarp2. **

	IF mapdelta GT 1 THEN BEGIN 
	   tempstrip(0,*) = stripn0(r,*)
	   tempstrip(1,*) = strp(r,*)
	   IF mapdelta GT (mapsize/360)*20 THEN BEGIN
		current = cmap(*,*,r)
		no_zero = WHERE (current GT 0)
		IF no_zero(0) LT 0 THEN fill = 1. ELSE $
			fill = median(current(no_zero))
		tempstrip(*,*) = fill		; ** Set 8/21/98 **
		print,'Fill =',fill
		printf,l1,'Fill =',fill
	   ENDIF
	   tempstrip(*,176:180) = max(cmap)>max(strp(r,*))
	   dummy = CONGRID(tempstrip,mapdelta*2,181,/interp)
	   newstrip = dummy(0:mapdelta,*)
	   IF k GT 0 THEN newstrip(0,*) = stripn0(r,*)
	   newstrip(mapdelta,*) = strp(r,*)
	
	   cmap(carrmap_x-mapdelta:carrmap_x,*,r) = newstrip(*,*)
	ENDIF ELSE cmap(carrmap_x,*,r) = strp(r,*) 
	
	stripn0(r,*) = strp(r,*)
	;if mapdelta gt 5 then stop
	 
;
;	** End carrwarp2.   **
	wset,0

	tvscl,cmap(*,*,r)<mapmax>mapmin,0,185*(r)
   ENDFOR
   IF mapdelta GT (mapsize/360)*20 THEN oldmedends=fltarr(20)
   oldmedends = shift(oldmedends,-1)
   IF carrmap_x GT 2 THEN oldmedends(19) = newmedend
   IF mapdelta GT (mapsize/360)*20 THEN oldstripmeds=fltarr(10)
   oldstripmeds = shift(oldstripmeds,-1)
   oldstripmeds(9) = newmed
   skiplast = 0

   skipped:
	   
;   IF k EQ numim-1 AND (t_end - t) GT 1/float(pixperday) THEN BEGIN
   IF k EQ numim-1 AND carrmap_x LT (mapsize-1) AND nextrot THEN BEGIN
   	FOR r=0,num_r-1 DO BEGIN               ; cycles through each radius
	   mapdelta = mapsize-1 - carrmap_x
	   IF mapdelta GT (mapsize/360)*20 THEN tempstrip(*,*)=0 ELSE BEGIN
	   	tempstrip(0,*) = strp(r,*)
	   	tempstrip(1,*) = strp(r,*)
	   ENDELSE
	   tempstrip(*,176:180) = max(cmap)
	   dummy = CONGRID(tempstrip,mapdelta*2,181,/interp)
	   newstrip = dummy(0:mapdelta,*)	   
	   newstrip(0,*) = strp(r,*)
	   cmap(carrmap_x:(mapsize-1),*,r) = newstrip(*,*)

	wset,0
	tvscl,cmap(*,*,r)<mapmax>mapmin,0,185*(r)
;stop
   	ENDFOR	; r loop
   ENDIF

   first = 0

ENDFOR     ;k loop

resultmaps = cmap
end_radius = radius
end_center(0) = cx
end_center(1) = cy
rollcor(1) = rotangle

free_lun,l1

IF nextrot THEN BEGIN
	carrmap_x = 0
	stripn0 = fltarr(num_r,181)
	oldstripmeds = fltarr(10)
	oldmedends = fltarr(20)
	stripends = fltarr(22)
	avg_strip_ends=[15,1,1.337]
	strp = fltarr(num_r,181)
	sta_center=fltarr(2)
	end_center=fltarr(2)
	cmap = fltarr(mapsize,181,num_r)
	tempstrip=fltarr(2,181)
	skiplast=0
	sta_radius=end_radius
	sta_center=end_center
ENDIF

save,filename=saveset, carrmap_x, y11,y22,y21,y12, x11,x22,x21,x12, strp, rad, stripn0, oldmedends, oldstripmeds, camnum, avg_strip_ends, cn, cmap, cam, mapsize, tempstrip, mapmax, mapmin, pixperday, skiplast, sta_radius, sta_center, rollcor, obs_timen0, mapdelta, platescl, cx,cy, radius

RETURN, reverse(resultmaps)

END ; carrmapmaker2