pro listmapmaker, listfile, num_r, radii, cmap, wlimb, display = display, range=range

;
;PURPOSE: 
;     Returns synoptic maps at multiple radii for a list of LASCO images
;
;INPUTS:
;	listfile   STR	name of file containing list of images -or-
;		   STRARR   array containing list of images
;	num_r	INT:	number of radii to make maps for
;	radii	FLTARR(num_r):		passed empty
;	cmap	FLTARR(mapsize,mapht,num_r):	passed empty
;	wlimb	STRING:	'wl' or 'el'; keywords for carrdate
;
;KEYWORDS:
;	DISPLAY	set to display images
;	RANGE	custom enter strip range (default is 180 to 0 degrees)
;
;
; OUTPUTS:	radii = fltarr(num_r)	:	radii of maps
;		cmap = fltarr(360,*,num_r):	carrington maps, one per radius
;		
;
;ROUTINES CALLED:
;    resize2.pro
;
; AUTHOR:	Nathan Rich, NRL, 1997/12/30
;
; MODIFIED:	NBR, 980203	make polar maps for K. Dere
;		NBR, 980210	call getc2c3norm2 (exposure time correction)
;		NBR, 990125	change platescl defn; cancel median check for PB
;				images
;

;COMMON warpit, carrmap_x,obs_timen0, stripn0, time0, duration, mapsize
COMMON map_disp, mapmax,mapmin, mapsize, cam, vert, cn, mapht, hdr, normalize, rat, startdate, platescl

sz = size (listfile)
numim=sz(1)
IF sz(0) EQ 0 THEN BEGIN
	filelist = readlist(listfile)
	sz = size(filelist)
	numim = sz(1)
ENDIF ELSE filelist = listfile
junk = lasco_readfits(filelist(0),hdr,/no_img)
startdate = hdr.date_obs+' '+hdr.time_obs
junk = lasco_readfits(filelist(numim-1),hdr,/no_img)
endate    = hdr.date_obs+' '+hdr.time_obs
cam = STRLOWCASE(hdr.detector)
;cn = 1922

IF keyword_set(RANGE) then BEGIN
	thetas = fltarr(2)
	print,'Enter start (top) and end (bottom) range of strips. '
	read, '(180 is South to West of image, 0 is North.): ',thetas
ENDIF ELSE thetas = [180.0,0.0]
sta_theta = thetas(0)
end_theta = thetas(1)
rad=radii

increm = (end_theta - sta_theta)/(mapht-1)
sz=size(radii)
num_r = sz(1)
skiplast= 0
factor=1.
;oldstripmeds = fltarr(10)
;oldmedends = fltarr(20)
;stripends = fltarr(22)
;avg_strip_ends=[15,3,2]
;skp=''
;near = ''
listfile = ''
endshort = 'no'
;val1=fltarr(num_r)
;val1(*) = 1
;val0= fltarr(num_r)
strip = fltarr(num_r,mapht)
;mnm = fltarr(1024,1024)		; c2 minimum
delta = 0.
first = 1
avg_med= [40,20,8]
carrmap_x = 0
stripn0 = fltarr(num_r,mapht)
obs_timen0 = STRARR(num_r)
tempstrip = FLTARR(2,mapht)
limb=fltarr(2)
IF strpos(wlimb,'wl') NE -1 THEN limb(0) = 1
IF strpos(wlimb,'el') NE -1 THEN limb(1) = 1

IF cn GT 1900 THEN BEGIN
   time0s = carrdate2(cn,WL=limb(0),EL=limb(1))
   time1s = carrdate2(cn+1,WL=limb(0),EL=limb(1))
ENDIF ELSE BEGIN
   ;print
   ;print,'Please enter start date of map. (First image is ',startdate,')'
   ;read,'(Ex.: 970414): ',sta_str
   ;read,'Please enter end date (0 hours):       ',end_str

   
   time0s = str2utc(startdate)
   time1s = str2utc(endate)
   cn  = utc2yymmdd(time0s)
ENDELSE

time0sec= utc2tai(time0s)
time1sec= utc2tai(time1s)
duration = (time1sec - time0sec)/86400.	; duration is in decimal days

time0_str = utc2str(time0s,/ECS)
time1_str = utc2str(time1s,/ECS)
t_start = time0s.mjd + time0s.time/86400000.
t_end = time1s.mjd + time1s.time/86400000.	

obs_timen0 = time0_str




FOR r=0,num_r-1 DO BEGIN
   REPEAT BEGIN
   IF cam eq 'c1' THEN BEGIN
	;READ,'Enter solar radius: ',radius
	;radii = [1.25,1.4,1.6,2]
	;radius = radii(r)
       rmin=1.05
       rmax=2
	;cx=255.2	;
	cx=254.2	;
	cy=247.7	;	These values are read from each header
	;cy=246		;		or computed from each header
	;rconv=87	;
	;platescl=11.2	;
	median0 = 0
	median1 = 0
	mapmax=1
	mapmin=-1
	disp_factor=1
	deltend=.8
	y11=390		;
	y12=410		;
	y21=100		;
	y22=120		;	These values are for the box normalization
	x11=85		;	    They are the vertices of two boxes
	x12=105		;
	x21=405		;
	x22=425		;
	camnum = 0
   ENDIF
   IF cam eq 'c2' THEN BEGIN
       READ, 'Enter solar radius between 2.2 and 5.2: ', radius
	;radii = [2.4,2.6,2.8,3,3.2,3.4,3.6,3.8,4,4.2,4.4,4.6,4.8,5,5.2,5.6,6]
	;radius = radii(r)
	rmin=2.2
	rmax=5.2
       cx=512
       cy=505
       ;rconv=79
	filt = 3
	platescl=12.1
	;median0 = 66
	;median1 = 85
	;median1= 2060
	;median0= 1700
	median0 = 1.03
	median1 = 1.1
	mapmax=4000
	mapmin=100
	IF rat EQ 'T' THEN BEGIN
		mapmax = 1.4
		mapmin =  .7
	ENDIF
	disp_factor=2
	;deltend=0.015
	deltend=0.05
	x11=10
	x12=70
	x21=960
	x22=1020
	IF t_start GT 50473 THEN BEGIN
		y11 = 850
		y12 = 890
		y21 = 130
		y22 = 170
	ENDIF ELSE BEGIN
		y11 = 758
		y12 = 798
		y21 = 225
		y22 = 265
	ENDELSE
	camnum = 1
   ENDIF
   IF cam eq 'c3' THEN BEGIN
      ;READ, 'Enter solar radius between 5 and 27: ', radius
      	radii = [5,7,10,15,20,24,28]
	radius = radii(r)
	rmin=5 
      rmax=30
      cx=518
      cy=530
      ;rconv=18
	filt=0
	platescl=56
	median0 = 1.3	; changed 8/6/98, NBR
	;median0 = 40	; changed 3/19/98, NBR
	median1 = 1.5	; changed 8/6/98, NBR
	mapmax=500	; changed 3/19/98, NBR
	mapmin=20
	;IF level EQ 'r' THEN BEGIN
	;	mapmax = 1.4
	;	mapmin = 1.2
	;ENDIF
	disp_factor=2
	deltend=.008	; changed 8/10/98
	;deltend = .025	; changed 8/26/98 for CR 1934
	IF cn LT 1930 THEN BEGIN
		x11=40
		x12=90
		x21=930
		x22=980
		y11 = 728
		y12 = 798
		y21 = 225
		y22 = 295
	ENDIF ELSE BEGIN
		x11=250
		x12=300
		x21=550
		x22=600
		y11 = 940
		y12 = 1010
		y21 = 10
		y22 = 80
	ENDELSE
	camnum = 2
   ENDIF
;   IF radius lt rmin OR radius gt rmax THEN $
;	PRINT, 'Radius is out of range. Please try again: '
   ENDREP UNTIL radius ge rmin and radius le rmax
   rad(r)=radius
ENDFOR
radii = rad		; for returning values
df = disp_factor

print
print,'time0 = ',time0_str,'    time1 = ',time1_str
print
print,'Found ',numim,'  images for CR ',cn, '   ',wlimb
print

pixperday = mapsize/duration		; factor to multiply days by to get 
					; the map size in x direction
					; (pixels per day)


OPENW,l1,cam+wlimb+'_'+TRIM(STRING(cn))+'_'+TRIM(STRING(rad(0)))+'.log' ,/get_lun
OPENW,l2,cam+wlimb+'_'+TRIM(STRING(cn))+'_'+TRIM(STRING(rad(0)))+'.info',/get_lun
printf,l1,'Found ',numim,' images for CR ',cn,'   ',wlimb
printf,l2,'White marks near the top of the image indicate where there is no data.'
printf,l2
printf,l2,'Found ',numim,' images for CR ',cn
printf,l2
printf,l1,'time0 = ',time0_str,'    time1 = ',time1_str
printf,l2,'carrmap_x is true column value: range is from carrmap_x to +mapdelta'
printf,l2,''


IF keyword_set(DISPLAY) THEN window,2,xsize=512,ysize=512
window,0,xsize=mapsize,ysize=mapht

FOR k = 0,numim-1 DO BEGIN
   wset,2
   print,'carrmapmaker: k = ',k
   toskip = 'n'
   filename = filelist(k)
   print, 'Reading file ',filename  
   nothing = lasco_readfits(filename,hdr,/no_img)
   utcdate = str2utc(hdr.date_obs)
   yymmdd = utc2yymmdd(utcdate)
   solar_ephem,yymmdd,radius=radius,/soho
   ;platescl = hdr.platescl		; no platescl values in PB images
   rconv = radius*3600/platescl
   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 or tn GT t_end  THEN BEGIN
	printf,l1,'Skipped ',filename,': Out of Order or too close to previous.'
	print,'Skipping ',filename,': Out of Order or too close to previous.'
	GOTO, skipped
   ENDIF

;
;  ** End carrwarp1.   **
;

;
; 	** Put your method of pre-processing your images here **
;	   such as 
;	im0 = LASCO_READFITS(filename,hdr)
;	   Remember to use 
;	im = REDUCE_STD_SIZE(im0,hdr,/full)
;	   as part of pre-processing.
;

   IF cam EQ 'c1' THEN BEGIN
  	IF level NE 'g' THEN im0 = getc1imwl(filename,obs_timen) $   ; ** white light C1
	ELSE im0 = getc1im3(filename,obs_timen)		; ** green line C1
	im = im0
	platescl = platescl/2
   ENDIF ELSE BEGIN
         IF (mapdelta EQ 0 AND carrmap_x NE 0) THEN BEGIN
	       print,'Skipping file ',filen
	    GOTO, skipped
        ENDIF
   	print, 'Reading file ',filename  
        IF rat EQ 'T' THEN im0 = getc2c3norm2(filename,hdr,/GET_BKG) $
	   ELSE im0 = getc2c3norm2(filename,hdr)
   	im = reduce_std_size(im0,hdr,/full)		; sizes im to 1024x1024 
	;platescl = platescl*(hdr.lebxsum > hdr.sumcol)
	;cx = cx/(hdr.lebxsum > hdr.sumcol)
	;cy = cy/(hdr.lebxsum > hdr.sumcol)
   ENDELSE
;
;	** END preprocess image **
;
   utcdate = str2utc(obs_timen)
   yymmdd = utc2yymmdd(utcdate)
   solar_ephem,yymmdd,radius=radius,/soho
   rconv = radius*3600/platescl
;   IF first THEN BEGIN
;	sta_radius = radius
;	sta_center(0) = cx
;	sta_center(1) = cy
;   ENDIF

   m = median(im)
   mx = max(im)

   IF keyword_set(DISPLAY) THEN BEGIN
	imsize= SIZE(im)

	IF imsize(1) GT 600 THEN BEGIN
	   xs = imsize(1)/2
	   ys = imsize(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
   	tvscl,hist_equal(disp_im)
;	IF k LT 3 THEN stop
   ENDIF

;  ** Begin box (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*m,null1)
   	ch2=where(squ2 NE 0 AND squ2 LT 1.4*m,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,3
		GOTO,skipped
   	ENDELSE		; factor is about im = 0.12*im for c2, 0.15*im for c3?
   ENDIF
	
;  ** End box normalization. **

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)	

wait,1

   IF m LE 0 and strpos(hdr.polar,'P') EQ -1 THEN BEGIN
	   	print,'Skipping file: Median out of range.'
	   	printf,l1,'Skipping file : Median out of range.'
	   	GOTO,skipped
   ENDIF


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

	IF vert EQ 'y' THEN BEGIN
	   xcoor = cx + westlimb*r0 + (1 + westlimb)*r0
	   strip(r,*) = congrid(im(xcoor,*),1,mapht)
	   IF keyword_set(DISPLAY) THEN disp_im(xcoor/disp_factor, *) =mx
	ENDIF ELSE BEGIN

	theta = sta_theta - 90
	FOR i=0,mapht-1 DO BEGIN
	   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=cy-r0*sin(th)
           strip(r,i)=INTERPOLATE(im,xcoor,ycoor)
	   IF keyword_set(DISPLAY) THEN BEGIN
	   	;disp_im(xecoor/disp_factor,ycoor/disp_factor)=mx
	   	;disp_im(xwcoor/disp_factor,ycoor/disp_factor)=mx
	   	disp_im(xcoor/disp_factor,ycoor/disp_factor)=mx
		;tvscl,disp_im
	   ENDIF
	   theta = theta+increm
	ENDFOR   ;	i loop
	ENDELSE

; 	** End take strip. **
;
;	w = where(stripn0)
;	IF w(0) LT 0 THEN stripn0(r,*) = strip(r,*)
;	;maxmin,strip(r,*)
;	strip(r,*) = strip(r,*)*avg_strip_ends(r)/newmedend
;	strip(r,*) = avg_med(r)*(strip(r,*)/median(strip(r,*)))
;
	IF median(strip(r,(mapht*0.2):(mapht*0.8))) LE 0 and strpos(hdr.polar,'P') EQ -1 $
	THEN BEGIN
		print,'Skipping file: strip out of range'
		GOTO, skipped
	ENDIF
   ENDFOR	; r loop


   IF keyword_set(DISPLAY) THEN BEGIN
	wset,2
   	tvscl,hist_equal(disp_im)
;	IF k LT 3 THEN stop
   ENDIF
	
   IF cam EQ 'none' 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.'
	   	GOTO,skipped
	   ENDIF
   	ENDELSE
   ENDIF

   carrmap_x = carrmap_x + mapdelta

   IF carrmap_x GT mapsize-1 THEN carrmap_x = mapsize -1
   IF mapdelta GT 1 THEN printf,l2,'carrmap_x =',359-carrmap_x,'    mapdelta =',mapdelta

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

	IF mapdelta GT 1 THEN BEGIN 
	   tempstrip(1,*) = strip(r,*)
	   tempstrip(0,*) = stripn0(r,*)
	   IF (deltaday GT 1.5 or total(stripn0(r,*)) LE 0) and strpos(hdr.polar,'P') EQ -1 $
	   THEN BEGIN
		current = cmap(*,*,r)
		no_zero = WHERE (current GT 0)
		IF no_zero(0) LT 0 THEN no_zero = WHERE (strip(r,*) GT 0)
		fill = median(current(no_zero))
		tempstrip(*,*) = fill		; ** Set 8/21/98 **
		print,'Fill =',fill
		printf,l1,'Fill =',fill
	   ENDIF
	   tempstrip(*,mapht-5:mapht-1) = max(cmap)>max(strip(r,*))
	   dummy = CONGRID(tempstrip,mapdelta*2,mapht,/interp)
	   newstrip = dummy(0:mapdelta,*)
	   IF k GT 0 and NOT(first) THEN newstrip(0,*) = stripn0(r,*)
	   newstrip(mapdelta,*) = strip(r,*)
	
	   cmap(carrmap_x-mapdelta:carrmap_x,*,r) = newstrip(*,*)
	ENDIF ELSE cmap(carrmap_x,*,r) = strip(r,*) 
	
	stripn0(r,*) = strip(r,*)
	;if mapdelta gt 5 then stop
	 
;
;	** End carrwarp2.   **
	wset,0

	tvscl,hist_equal(cmap(*,*,r)),0,(mapht+4)*(r)
   ENDFOR

   IF k EQ numim-1 AND carrmap_x LT (mapsize-1) 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,*) = strip(r,*)
	   	tempstrip(1,*) = strip(r,*)
	   ENDELSE
	   tempstrip(*,176:180) = max(cmap)
	   dummy = CONGRID(tempstrip,mapdelta*2,mapht,/interp)
	   newstrip = dummy(0:mapdelta,*)	   
	   newstrip(0,*) = strip(r,*)
	   cmap(carrmap_x:(mapsize-1),*,r) = newstrip(*,*)

	wset,0
	tvscl,hist_equal(cmap(*,*,r)),0,(mapht+4)*(r)
;stop
   	ENDFOR	; r loop
   ENDIF


   obs_timen0 = obs_timen
   first = 0

   skipped:
ENDFOR     ;k loop

free_lun,l1
free_lun,l2
cmap = reverse(cmap)
print
print,'time0 = ',time0_str,'    time1 = ',time1_str
print

END ; listmapmaker