;+ ;Procedure: CARRMAPMAKER2 ; $Id: carrmapmaker2.pro,v 1.15 2017/12/19 23:01:47 nathan Exp $ ; ;PURPOSE: ; Returns carrington maps at multiple radii for a passed ; carrington number at radii read from user in this routine. ; Program only works for maps beginning after ; January,6 1996. ; ;INPUTS: ; cn INT: carrington number ; cam STRING: camera, i.e. 'c2' ; num_r INT: number of radii to make maps for ; rad FLTARR(num_r): radii for maps ; cmap FLTARR(mapsize,ht,num_r): passed empty ; pol STRING: 'p','u' or ''; indicates polarizer value ; limb INTARR: 0 or 1 tells which limb for carrdate procedure ; level STRING: 'ql' for quick-look, 'r' for doing ratio ; disp INT: 1 to display images ; wlimb STRING: 'wl' or 'el' ; hdr STRUCT: dummy header sent back to carrmap3 for fits header ; ;KEYWORDS: ; FULL360 If set, go all the way around, not just 180 deg. ; MAP_SIZE IF set, then use it's value for horizontal size ; DEBUG= Number of seconds to wait before starting next image ; This routine can ask the user for each radius. ; ; OUTPUTS: rad = fltarr(num_r) : radii of maps ; cmap = fltarr([360 or 720],ht,num_r): carrington maps, one per radius ; ; ;ROUTINES CALLED: ; qdb.pro ; carrdate2.pro ; resize2.pro ; getc2c3norm2.pro ; ; AUTHOR: Nathan Rich, NRL, Nov. 1996 ; Julia Kraemer, NRL, June 7, 1996 ; ; MODIFIED: ; $Log: carrmapmaker2.pro,v $ ; Revision 1.15 2017/12/19 23:01:47 nathan ; commit current versions ; ; Revision 1.14 2015/05/15 18:47:38 nathan ; need /pivot in rot() call ; ; Revision 1.13 2015/03/26 15:23:40 nathan ; misc improvements ; ; Revision 1.12 2013/03/29 20:20:14 nathan ; bump up ss (strip to use for test) ; ; Revision 1.11 2013/03/29 19:04:41 nathan ; call mk_img with /PAN so all comes out as 1024; use fhdr for hdr ; ; Revision 1.10 2011/09/09 20:30:39 nathan ; extend common get_im ; ; Revision 1.9 2010/12/09 22:42:40 nathan ; adjustments Jul 21 2010 ; ; Revision 1.8 2010/05/18 15:58:04 nathan ; fix off_hdr bug ; ; Revision 1.7 2010/05/18 15:56:31 nathan ; fix off_hdr bug ; ; Revision 1.6 2010/05/18 15:40:14 nathan ; change roll and toskip logic ; ; Revision 1.5 2010/02/11 16:11:03 nathan ; use absolute strip median cutoff value for c1 ; ; Revision 1.4 2010/02/08 22:34:19 nathan ; Numerous mods for C1 maps including: use NaN not white to indicate missing data; adjustments to ; skip criteria based on image and strip values; enhanced logging to l1; use 9-pixel median instead ; of interpolate for strip values; added bell for skipped images ; ; 961129 by N. RIch each image subtends one pixel; checks for other ; values of val0 if 0 when k=0 ; 961203 by N. RIch if more than one pixel between images, use ; congrid to fill space ; 961206 by N. Rich Create log file; rotate images from 960521; ; Check median = 85 (for 1024x576 images) ; 961209 by N. Rich change CONGRID statement; moved carrmap_x ; increment statement to after validity check ; 961210 by N. Rich print time0 and time1; if mapdelta GT 10, don't ; do congrid ; 961211 by N. Rich skip bad files ; 961213 by N. Rich add level variable; if mapdelta GT 20, don't do ; congrid; check for forward time ; 961217 by N. Rich add disp variable for displaying image ; 970102 by N. Rich make time units consistent ; 970103 by N. Rich rotate images for Nov. 21,22 ; 970103 by N. Rich add white strip to top of gaps ; ; 970108 by J. Kraemer add east limb option ; 970123 by N. Rich changed to val0(r)=median(strip(115:135)) ; 970213 by N. Rich changed check for Out of Order ; 970226 by N. Rich do normalization for whole image rather than ; each strip ; 970311 by N. Rich make map display scale vary with camera ; 970402 by N. RIch add empty space to final columns if no images ; 970414 by N. Rich accept c1 images for Clarence ; 9705 by N. Rich use c1 images generally ; 970714 by N. Rich change criteria for data gaps ; 970715 by N. Rich fix final gap coverer ; 970722 by N. Rich check each strip against previous and reject ; image if it is too different ; 970730 by N. Rich sort s by date_obs ; 970801 by N. Rich skip image if mapdelta LE 0 ; 970807 by N. Rich scan strip ends; move imdisplay part around ; 970827 by N. Rich change scale; modify name of log files ; 970915 by N. Rich add get_im common block, compute ff_ratio; change ; cx,cy for c1 ; 971021 by N. Rich compute rconv for each image separately for c2 and c3 ; 971103 by N. Rich prompt for mapsize ; 971107 by N. Rich move log print line; prompt for median skip; don't ; check image median ; 971110 by N. Rich add skiplast variable ; 971114 by N. Rich use strip files if available ; 971216 by N. Rich make automatic (auto-enter radii, mapsize) ; 971217 by N. Rich if error in readfits, enter in log ; 980212 by N. Rich added variable boxes for normalization ; 980312 by N. Rich cancel box normalization; auto choose mapsize ; 980313 by N. Rich do not multiply image by constant (of about 0.12) ; 980413 by N. Rich use m = median(im(*,y21:y12)) ; 980514 by N. Rich change upper limit of strip median, ends ; 980610 by N. Rich extend warpit COMMON block for fits header ; 980612 by N. Rich pass rad (radii) already full ; 980629 by N. Rich use reduce_std_size instead of resize2 ; 981216 by N. Rich Use ht variable instead of 181 for mapheight ; 990114 by N. Rich Add FULL360 keyword ; Feb 2000 by N. RIch Remove source from query; change GET_BKG init.; ; Base level on LZ keyword, levelstr; Replace strip with ; strp to avoid IDL conflict ; Apr 2000, NBR, Add half a column's time to t_start; always use first image (k=0) ; May 2000, NBR, Use mk_img instead of getc123; don't use avg_strip_ends ; Aug 2000, NBR, Add MAP_SIZE keyword and take it out of common block; compare strip,stripends ; 2001.11.05, NBR - Add SCCS version tag ; 2002.03.14, NBR - Add skipwait, medianhist, variable sunc from mk_img, roll in mk_img ; 2002.09.03, nbr - Expand logging, add l1 to COMMON get_im ; 04.05.05, nbr - fix C1 center, tinker with disp_factor and df and ffv ; 06.04.25, nbr - Add filetype so does not pick up Level-1 images; remove date_mod ; from query for LZ images ; 06.10.03, nbr - Make missing-bar min. of 5 ; 2008.03.27, nbr - Accomodate new qdb output structure ; 09.07.07, nbr - Changed datatype of normalize ; 2014.12.04, nbr - Add /NORMALIZE option; always wait before skipping image ; ; %W% %H% - IDL NRL LASCO Library ; ;- ; pro carrmapmaker2, cn, cam, num_r,rad, cmap,pol,limb,level,disp, wlimb, dhdr, FULL360=full360, RUNNING_DIFF=running_diff, LZ=lz, MAP_SIZE=map_size, DEBUG=debug, $ NORMALIZE=normalize COMMON warpit, carrmap_x, obs_timen0, stripn0, time0, duration, rat, $ time0_str, time1_str, box_ref, sta_radius, end_radius, $ sta_center, end_center, rollcor, platescl, ht, addfac, numim, $ ctr, processhdr, history1, maplog, avtrend COMMON map_display, mapmax,mapmin, debugtime0 COMMON get_im, off_s,off_times,ff_ratio, off_hdr, prev, pylonim, l1, ffv,ref_exp IF keyword_set(DEBUG) THEN nomwait=debug ELSE nomwait=0 debug=0 nomroll=0. maxmissing=0.11 ; fraction IF datatype(box_ref) EQ 'UND' THEN box_ref=0 IF datatype(debugtime0) EQ 'UND' THEN debugtime0=300 addfac=0. fill = 1.0 ; ** set 6/20/00 barht = ht/36 < 20 prev = fltarr(1024,1024) rotangle = 0 skiplast= 0 skipwait= 1 factor=1. medianhist = fltarr(10) oldstripmeds = fltarr(10) oldmedends = fltarr(20) stripends = fltarr(22) avg_strip_ends=[15,1,1.337] skp='' near = '' listfile = '' sta_str = '' end_str = '' endshort = 'no' IF cam EQ 'c1' THEN scale = 60 ELSE scale = 200 mscale = 20 val1=fltarr(num_r) val1(*) = 1 val0= fltarr(num_r) strp = fltarr(num_r,ht) mnm = fltarr(1024,1024) ; c2 minimum delta = 0. first = 1 sta_center=fltarr(2) end_center=fltarr(2) logroot=cam+wlimb+'_'+TRIM(STRING(cn))+'_'+strmid(TRIM(STRING(rad[0])),0,3) OPENW,l1, logroot+'.log' ,/get_lun OPENW,l2, logroot+'.info',/get_lun printf,l1,'Times are local time: systime().' IF level EQ 'gm' THEN narrow=0 ELSE narrow=1 ; whether to include next lower wavelength IF cn GT 1900 THEN BEGIN 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)) ENDIF ELSE BEGIN print, 'Carrington number invalid. Please start over.' return ENDELSE time0sec= utc2tai(time0s) time1sec= utc2tai(time1s) duration = (time1sec - time0sec)/86400. ; duration is in decimal days IF duration GT 28 THEN BEGIN help,duration read,'Return to continue>',skp ENDIF time0_str = utc2str(time0s,/ECS) time1_str = utc2str(time1s,/ECS) t_start = time0s.mjd + double(time0s.time/86400000.) t_end = time1s.mjd + double(time1s.time/86400000.) obs_timen0 = time0_str ;IF cn LT 1921 THEN ss=1 ELSE ss=num_r-2 ss=(num_r/2)>2 ;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 disp_factor=0.5 ; display is fraction of images used ffv=1 ; uses binned images ;cx=255.2 ; original values ;cy=246 ; ;cx=254.2 ; header values ;cy=247.7 ; ;cx=511./(2/(ffv+1)) ; fudged to center occulter, 11/23/98 ;cy=494./(2/(ffv+1)) ; ;platescl=5.8 ; fudged to put occulter at 1.07 R, 11/23/98 ;rconv=87 ;platescl=11.74 ; computed by NBR, 11/19/98 platescl=5.8 * (2/(ffv+1)) ; defined platescl, 04.05.06, nbr mapmax=7 mapmin=-1 deltend0=.8 y11=720/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y12=820/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y21=200/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y22=300/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x11=210/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x12=310/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x21=710/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x22=810/(2/(ffv+1)) ; FFV=1 means 1024x1024 images camnum = 0 addfac=3.0 ; since FeXIV C1 is a difference image, result is biased negative stripmin=addfac+1. stripmax=addfac*6+1 median0 = addfac ; after addfac and before normalization median1 = median0 ENDIF IF cam eq 'c2' THEN BEGIN nomroll=0.50 ;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) disp_factor=0.5 ; display is half size of images used ffv=1 ; uses 1024x1024 images rmin=2.2 rmax=5.2 cx=512./(2/(ffv+1)) cy=505./(2/(ffv+1)) ;rconv=79 filt = 3 platescl=12.1 * (2/(ffv+1)) ;median0 = 66 ;median1 = 85 ;median1= 2060 ;median0= 1700 median0 = 1.03 ; This value is for screening first image only median1 = 1.58 ; cosmic ray storm case mapmax=4000 mapmin=100 IF rat NE 'F' THEN BEGIN mapmax = 1.8 mapmin = .8 ENDIF ;deltend=0.015 ;deltend=0.05 deltend0=0.22 ;2/2/00 x11=200/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x12=260/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x21=600/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x22=660/(2/(ffv+1)) ; FFV=1 means 1024x1024 images IF t_start GT 51910 THEN BEGIN ; 2001/01/01 y11 = 960/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y12 = 1000/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y21 = 20/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y22 = 60/(2/(ffv+1)) ; FFV=1 means 1024x1024 images ENDIF ELSE BEGIN y11 = 758/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y12 = 798/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y21 = 225/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y22 = 265/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x11=20/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x12=80/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x21=960/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x22=1020/(2/(ffv+1)) ; FFV=1 means 1024x1024 images ENDELSE camnum = 1 med_expt = 25 ENDIF IF cam eq 'c3' THEN BEGIN nomroll=-0.23 ;READ, 'Enter solar radius between 5 and 27: ', radius ;radii = [5,7,10,15,20,25] ;radius = radii(r) disp_factor=0.5 ; display is half size of images used ffv=1 ; uses 1024x1024 images rmin=5 rmax=27 cx=518./(2/(ffv+1)) ; FFV=1 means 1024x1024 images cy=530./(2/(ffv+1)) ; FFV=1 means 1024x1024 images ;rconv=18 filt=0 platescl=56*(2/(ffv+1)) ; FFV=1 means 1024x1024 images median0 = 1.0 ; changed 7/6/99, NBR ;median0 = 40 ; changed 3/19/98, NBR median1 = 1.0 ; changed 7/6/99, NBR mapmax=1.1 ; changed 7/6/99, NBR mapmin=.95 IF rat NE 'F' THEN BEGIN mapmax = 1.15 mapmin = 0.9 ENDIF ;deltend=.015 ; changed 2/8/00, CR 1949 EL ;deltend = .025 ; changed 8/26/98 for CR 1934 deltend0 = .05 ; CR 1905 IF cn LT 1930 THEN BEGIN x11=40/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x12=90/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x21=930/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x22=980/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y11 = 623/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y12 = 693/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y21 = 260/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y22 = 330/(2/(ffv+1)) ; FFV=1 means 1024x1024 images ENDIF ELSE BEGIN x11=250/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x12=300/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x21=550/(2/(ffv+1)) ; FFV=1 means 1024x1024 images x22=600/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y11 = 940/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y12 = 1000/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y21 = 60/(2/(ffv+1)) ; FFV=1 means 1024x1024 images y22 = 120/(2/(ffv+1)) ; FFV=1 means 1024x1024 images ENDELSE camnum = 2 med_expt = 20 ss=3 ENDIF IF keyword_set(RUNNING_DIFF) THEN BEGIN mapmax=20 mapmin=0 median0 = 0.9 median1 = 1.1 deltend0 = 0.5 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 ;rad = radii df = disp_factor avg_med= [40,20,8] carrmap_x = 0 stripn0 = fltarr(num_r,ht) obs_timen0 = STRARR(num_r) tempstrip = FLTARR(2,ht) ;IF level EQ 'ql' THEN BEGIN ; levelnum = 1 ; levelstr = 'QUICK-LOOK' ;ENDIF ELSE BEGIN ; levelnum = 2 ; levelstr = 'LZ' ;ENDELSE IF keyword_set(LZ) THEN BEGIN levelnum='=2' qldupchk='' ENDIF ELSE BEGIN levelnum='=1' ;qldupchk='img_files.date_mod=img_leb_hdr.date_mod and ' qldupchk='' ENDELSE print print,time0_str print,'to' print,time1_str now=systime() IF cam EQ 'c1' THEN BEGIN polnum=0 IF narrow THEN BEGIN on_wl = 5302.4 width = 0.1 ENDIF ELSE BEGIN on_wl = 5302.0 width = 0.5 ENDELSE print print,'Using ',on_wl,' A for on-line.' print on_query = 'select img_files.filename, date_obs, fp_wl_upl, lp_num, naxis1, naxis2, diskpath from img_leb_hdr, img_files where img_leb_hdr.fileorig=img_files.fileorig and '+qldupchk+'camera=' +TRIM(STRING(camnum))+' and polar=' +TRIM(STRING(polnum))+' and filetype=1 and source=2 and (lp_num=8 or lp_num=9) and date_obs between "'+time0_str+'" and "'+time1_str+'" and fp_wl_upl between ' +TRIM(on_wl)+' and ' +TRIM(on_wl+width)+' and exptime < 50 and naxis1 > 300 and naxis2 > 256' s=qdb('lasco',on_query) print,'on_query=',on_query print printf,l1,'on_query=',on_query printf,l1 ;off_sz = size(off_s.date_obs) on_sz = size(s.date_obs) num_on = on_sz(1) numim=num_on ;numim = trim(string(num_off))+' offline, '+trim(string(num_on))+' online' goto, skipoff off_times=dblarr(num_off) IF level EQ 'g' THEN BEGIN FOR kk=0,num_off-1 DO BEGIN off_utc=str2utc(off_s.date_obs(kk)) off_times(kk)=utc2tai(off_utc) ENDFOR ENDIF lz05 = '/net/corona/cplex2/lz_data/level_05/' ffonim = lasco_readfits(lz05+'960910/c1/12014004.fts',ffon_hdr) ffoffim =lasco_readfits(lz05+'960910/c1/12014005.fts',ffoff_hdr) ffonim = REDUCE_STD_SIZE(TEMPORARY(ffonim), ffon_hdr,/full) ;** make all images 1024x1024 ffoffim = REDUCE_STD_SIZE(TEMPORARY(ffoffim), ffoff_hdr,/full) ;** make all images 1024x1024 ffoffim=ffoffim+1 ff_ratio = rebin( float(ffonim)/float(ffoffim), 1024/(2/(ffv+1)), 1024/(2/(ffv+1)) ) processhdr="online - offline*[960910/c1/(12014004.fts/12014005.fts)]" skipoff: ENDIF ELSE BEGIN IF (pol EQ 'p' OR pol EQ 'u') AND cam EQ 'c2' THEN polnum = 2 ELSE $ IF (pol EQ 'p' OR pol EQ 'u') AND cam EQ 'c3' THEN polnum = 3 ELSE $ polnum = 0 query = 'select img_files.filename, img_leb_hdr.date_obs, img_leb_hdr.exptime, img_leb_hdr.fp_wl_upl, img_leb_hdr.lp_num, img_leb_hdr.naxis1, img_leb_hdr.naxis2, img_files.diskpath from img_leb_hdr, img_files where img_leb_hdr.lp_num!=17 and img_leb_hdr.lp_num!=5 and img_leb_hdr.lp_num!=7 and '+qldupchk+' img_leb_hdr.fileorig=img_files.fileorig and img_leb_hdr.camera=' +TRIM(STRING(camnum))+' and img_files.filetype=1 and img_files.source ' +levelnum+' and img_leb_hdr.polar=' +TRIM(STRING(polnum))+' and img_leb_hdr.filter=' +TRIM(STRING(filt))+' and img_leb_hdr.naxis1>511 and img_leb_hdr.naxis2>511 and img_leb_hdr.date_obs between "' +time0_str+'" and "'+time1_str+'" and img_leb_hdr.exptime between ' +trim(string(0.75*med_expt))+' and '+trim(string(1.25*med_expt)) ; 11/13/14, nr - removed "and img_leb_hdr.naxis1=img_leb_hdr.naxis2 and (img_leb_hdr.lp_num=9 or img_leb_hdr.lp_num=14) " from query because skipping c2 files in late 1998 qu=query t0=systime(1) ;seconds s=qdb('lasco',qu) t1=systime(1) print,query help,/struct,s wait,3 IF (datatype(s) NE "STC") THEN BEGIN printf,l1,"cn ",cn," No Data Found" printf,l2,"cn ",cn," No Data Found" wait,10 free_lun,l1 free_lun,l2 return ENDIF sz = SIZE(s.filename) numim = sz(1) printf,l1,'Query = '+qu printf,l1 printf,l1,now+': query took '+STRING(t1-t0)+' seconds.' ENDELSE print print,'time0 = ',time0_str,' time1 = ',time1_str print print,'Found ',numim,' images for CR ',cn, ' ',wlimb print printf,l1,'Found ',numim,' images for CR ',cn,' ',wlimb printf,l1,'time0 = ',time0_str,' time1 = ',time1_str printf,l1 ;print,'Sorting s.' ord = sort(s.date_obs) ;read,'Enter mapsize (720 or 360): ',mapsize ta=str2utc(s(ord(0)).date_obs) tb=str2utc(s(ord(numim-1)).date_obs) ta2b=float(tb.mjd-ta.mjd) files=s(ord).diskpath+'/'+s(ord).filename obs_times = s[ord].date_obs IF keyword_set(MAP_SIZE) THEN mapsize=map_size ELSE $ IF cam EQ 'c2' or cam EQ 'c3' THEN IF numim/ta2b GT 28 THEN mapsize = 720 ELSE mapsize = 360 help,mapsize map_size=mapsize cmap=FLTARR(mapsize,ht,num_r) pixperday = mapsize/duration ; factor to multiply days by to get ; the map size in x direction ; (pixels per day) printf,l2,'NaN strips 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,l2,'carrmap_x is true column value: range is from carrmap_x to +mapdelta' printf,l2,'' IF disp EQ 1 THEN window,2,xsize=512,ysize=512 window,0,xsize=mapsize,ysize=ht+10,title='CR '+trim(cn)+' '+wlimb+' R='+trim(rad[ss]) ctr=0 ;stop FOR k = 0,numim-1 DO BEGIN print,'carrmapmaker: k = ',k nooverride=1 debugtime=debugtime0 toskip = 'n' rotangle = 1 ; initialize for call to mk_img sunc = 1 filename = s(ord(k)).filename source = strmid(filename,1,1) filen = files[k] obs_timen = obs_times[k] IF filename eq '22001176.fts' or $ filename eq '22001193.fts' or $ filename eq '22001196.fts' or $ filename eq '22007037.fts' or $ filename eq '22009079.fts' or $ filename eq '22013010.fts' or $ filename eq '22009417.fts' or $ filename eq '22037273.fts' or $ filename eq '22529923.fts' or $ filename eq '32004565.fts' or $ filename eq '32004570.fts' or $ filename eq '32004574.fts' or $ filename eq '32004575.fts' or $ filename eq '32004578.fts' or $ filename eq '32014488.fts' or $ filename eq '32004579.fts' THEN GOTO, skipped IF filename EQ '22103331.fts' THEN stop ; ** Begin carrwarp1. ** ; dayn = str2utc(obs_timen) tn = dayn.mjd + DOUBLE( dayn.time/86400000.) ; time is milliseconds, changed to day units ; day is modified julian day tn0 = t_start + (0.5/pixperday) + carrmap_x/pixperday deltaday = tn - tn0 ; time between last column and present image fltmapdelta=deltaday*pixperday mapdelta = ROUND(deltaday*pixperday) print, 'tn =',tn,' deltaday =',deltaday,' mapdelta =',mapdelta ;stop IF mapdelta LE 0 and NOT(first) THEN BEGIN printf,l1,strmid(obs_timen,0,10)+' Skipped '+filename+': mapdelta= '+trim(string(mapdelta)) print,'Skipping ',filename,': Out of Order or too close to previous.' GOTO, skipped ENDIF ; ; ** End carrwarp1. ** skipwait=skipwait-1 number = strmid(filename,0,8) utcday = str2utc(obs_timen) yymmdd = utc2yymmdd(utcday) stripfile = findfile('/net/corona/export/raid2/ew_strips/' +yymmdd +'/' +cam +'/' +number +'_' +wlimb +'.fts') print,filen IF stripfile(0) EQ 'd' THEN BEGIN print,'Using stripfile...' printf,l1,obs_timen+' '+filename+': Using stripfile: ',stripfile strp= readfits(stripfile(0),hdr) ENDIF ELSE BEGIN ;*************************************************************** IF cam NE 'c1' THEN BEGIN ; IF (mapdelta EQ 0 AND carrmap_x NE 0) THEN BEGIN ; print,'Skipping file ',filen ; GOTO, skipped ; ENDIF ; IF pol EQ 'p' THEN im = getc2c3pol2(filen) ELSE $ ; IF pol EQ 'u' THEN im = getc2c3unpol2(filen) ELSE $ print, 'Reading file ',filen IF rat EQ 'ANY_YEAR' THEN get_bkg = 1 ELSE $ IF rat EQ 'T' THEN get_bkg = 2 ELSE get_bkg = 0 IF rat EQ 'F' THEN dorat=0 ELSE dorat=1 REPEAT BEGIN ;im0 = getc123(filen,hdr,GET_BKG=get_bkg, RUNNING_DIFF=running_diff) im0=mk_img(filen,0.8,1.3,fhdr,/automax, /NO_DISP, RATIO=dorat,USE_MOD=get_bkg, FIXG=0, $ RUNNING=running_diff, ROLL=0, SUNC=sunc, PAN=1, EXPFAC=efac) ;+++ mk_img: subtract bias; divide by exposure duration; divided by monthly model; fix binning IF im0[0] EQ -1 THEN BEGIN print,'READFITS ERROR? -- stopped. Type .cont to continue with next file.' stop im0=intarr(1024,1024) sunc={xcen:0,ycen:0} ENDIF ELSE $ hdr=lasco_fitshdr2struct(fhdr) ENDREP UNTIL im0(0) NE -1 im = im0 cx = sunc.xcen cy = sunc.ycen offinfo='' ENDIF ELSE BEGIN ;filen = filelist(k) IF strmid(level,0,1) NE 'g' THEN BEGIN im = getc1imwl(filen,obs_timen) ; ** white light C1; hdr in common block processhdr='getc1imwl.pro' ENDIF ELSE BEGIN im0 = make_image_c1(filen,hdr, OFF_HDR=off_hdr) ; ** FeXIV green line C1 processhdr='fpc1_2img with fVof0996.fts and fVon0996.fts' ENDELSE cx=hdr.crpix1 cy=hdr.crpix2 fp_wl_upl=hdr.fp_wl_up im = (im0+addfac) offinfo=', off: '+off_hdr.filename+' @ '+off_hdr.date_obs+' '+strmid(off_hdr.time_obs,0,8) ENDELSE lpnum=cnvrt_lp(hdr.lp_num) 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 szim = SIZE(im) IF szim(0) EQ 0 THEN BEGIN print,'File not found--skipping' printf,l1,obs_timen+' '+filename+'File not found--skipping. time='+systime() GOTO, skipped ENDIF IF szim(1) GT 600 THEN m = median(im(*,y21:y12)) $ ELSE m=median(im(*,y21/2:y12/2)) mx = max(im) IF skipwait GE 0 or $ (tn LT 51746 and tn GT 51739) or $ ; 14 Jul 00 event (tn LT 51861 and tn GT 51856) or $ ; 9 Nov 00 event (tn LE 52220 and tn GE 52217) $ ; 4 Nov 01 Event THEN skipflag='n' ELSE skipflag='y' medianok=1 IF noT(keyword_set(RUNNING_DIFF)) and $ ; s(ord(k)).naxis2 GT 512 and $ skipflag NE 'n' THEN BEGIN IF (m GT 1.1*median0 or m LT 0.9*median0) THEN BEGIN $ toskip='y' ; read,'Skip this one? [y/n]',skip print,'Median = ',m print,medianhist print,string(7b) ; bell print,'Skipping file: Median out of range.' medianok=0 debugtime=30 nooverride=1 ENDIF ENDIF rotangle=hdr.crota crota=-hdr.crota+nomroll IF abs(hdr.crota-nomroll) GT 1. THEN BEGIN ; SC roll of -2.1 on 8/5/96 19:30 is apparently an error print,string(7b) ; bell message,'Rotating image '+trim(crota)+' deg',/info printf,l1,obs_timen+' '+filename+' Rotating image '+trim(crota)+' deg' wait,2 im=rot(temporary(im),crota,1.,cx,cy,/pivot) ENDIF ;else wait,nomwait IF (tn GT 50224.0 AND tn LT 50224.83) OR $ (tn GT 50408.724 AND tn LT 50409.33) then BEGIN ;im = ROTATE(im,3) ; rotate images from 960521 and 961122 IF cam EQ 'c3' THEN im = SHIFT(im,-12,25) ENDIF ;IF tn LT 51564 THEN rotangle = -1*180*GET_SOLAR_ROLL(hdr,/medi)/!pi $ ;ELSE BEGIN ; pnt = GET_SC_POINT(obs_timen) ; rotangle = -1*pnt.sc_roll ; degrees ;ENDELSE ; ** Now getting rotangle from mk_img.pro ; ;help,rotangle ; ** roll done in mk_img ; IF ABS(rotangle) GT 1 THEN BEGIN ; im =rot(temporary(im),rotangle,1,cx,cy,/pivot,/interp) IF carrmap_x LE 1 THEN $ IF cam eq 'c1' THEN $ rollcor[0]=get_roll_or_xy(hdr,'roll',/deg) ELSE $ rollcor(0)=rotangle ; print,'Rotating ',filename,rotangle,' degrees CW.' ; printf,l1,obs_timen+' '+filename+': Rotating ',rotangle,' degrees CW.' ; ENDIF ; IF k eq 1 AND cam eq 'c2' THEN BEGIN ; mnm=lasco_readfits('/net/corona/cplex3/monthly/c2min.fts',hm) ; bias=470 ; exp=hm.EXPTIME ; mnm=(im-bias)/exp ; ENDIF IF disp EQ 1 THEN BEGIN imsize= SIZE(im) IF imsize(1) GT 600 THEN BEGIN xs = imsize(1)*disp_factor ys = imsize(2)*disp_factor disp_im = rebin(im,xs,ys) ENDIF ELSE disp_im = im ;IF cam EQ 'c1' or cam EQ 'c3' THEN BEGIN ;IF normalize GE 0 THEN BEGIN 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 ;ENDIF wset,2 IF pol EQ 'p' OR pol EQ 'u' THEN scale = 80 tvscl,disp_im>mapmin1),null1) ch2=where(squ2 GT 0 AND abs(squ2) LT 5*(addfac>1),null2) percentbad = 1 - float(null1+null2)/nsqus IF percentbad GT 0.5 THEN BEGIN print,string(7b) ; bell print,'Skipping file: ',fix(100*percentbad),'% of boxes out of range.' printf,l1, obs_timen+' '+ filename+' '+ trim(fix(100.*percentbad))+ '% of boxes out of range at Carrmap_x='+string(carrmap_x,'(i3)')+ ', LP='+string(lpnum,'(i2)')+ offinfo wait,2 GOTO,skipped ENDIF ;IF cam EQ 'c1' THEN BEGIN ; add factor to get rid of negative values ; squ1=(squ1+addfac)>0 ; squ2=(squ2+addfac)>0 ;ENDIF IF float(null1)/nsqu LT 0.6 THEN a1=0 IF float(null2)/nsqu LT 0.6 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 mdn1 = median(squ1) mdn2 = median(squ2) int_corr=a1*mdn1 + a2*mdn2 factor = box_ref/int_corr ;** box_ref 0.17 for c1? print,'factor = ',factor,' toplft =',null1,' botrt =',null2 ;if carrmap_x gt 218 then stop IF int_corr GT 0 THEN BEGIN IF keyword_set(NORMALIZE) THEN im=factor*im ELSE message,'Not normalizing images.',/info ENDIF ELSE BEGIN print,string(7b) ; bell print,'Skipping file: box median out of range.' printf,l1,obs_timen+' '+filename+' 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') IF (westlimb) and carrmap_x eq 406 and cn eq 2164 then stop IF keyword_set(FULL360) THEN range = 360. ELSE range = 180. 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,ht-1 DO BEGIN theta=i*(range/(ht-1))-90.0 th=theta*!PI/180. IF westlimb NE -1 or keyword_set(FULL360) 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 IF xwcoor GE 0 and xwcoor LE 1023 and $ xecoor GE 0 and xecoor LE 1023 and $ ycoor GE 0 and ycoor LE 1023 THEN BEGIN ;if r eq 0 then print,xcoor,ycoor ;strp(r,i)=INTERPOLATE(im,xcoor,ycoor) strp[r,i]=median(im[xcoor-1:xcoor+1,ycoor-1:ycoor+1]) samplemethod='Sample method: 9-pixel median' ;if i eq 100 and r eq 2 then stop IF disp EQ 1 THEN BEGIN ydisp=ycoor*disp_factor xdisp=xcoor*disp_factor xedisp=xecoor*disp_factor xwdisp=xwcoor*disp_factor ;disp_im(xedisp,ydisp)=mx disp_im(xdisp,ydisp)=mx ENDIF ENDIF ; IF k eq 1 THEN mstrip(i)=interpolate(mnm,xcoor,ycoor) ENDFOR ; i loop ; ** End take strip. ** ; ; w = where(stripn0) ; IF w(0) LT 0 THEN stripn0(r,*) = strp(r,*) ; ;maxmin,strp(r,*) ; strp(r,*) = strp(r,*)*avg_strip_ends(r)/newmedend ; strp(r,*) = avg_med(r)*(strp(r,*)/median(strp(r,*))) ; ENDFOR ; r loop ; IF source EQ 2 and cam NE 'c1' THEN savestrips,strp,hdr,factor,rad,wlimb ;IF k LT 2 THEN stop print,obs_timen+' '+filename+': carrmap_x ='+STRING(carrmap_x,'(i3)')+', rconv='+STRING(rconv,'(f7.2)') +' median(im) ='+STRING(m,'(f6.2)')+' max(im) =' +STRING(mx,'(f8.2)') help,lpnum IF disp EQ 1 THEN BEGIN wset,2 IF pol EQ 'p' OR pol EQ 'u' THEN scale = 80 tv,bytscl(disp_im,min=mapmin,max=mapmax) xyouts,10,10,obs_timen,size=2,/device ENDIF ;IF k LT 54 and k GT 44 THEN stop ;****************************************************************** ENDELSE print ;; ** BEgin check strip ** sst = SORT(strp[2,*]) lowstrp = strp[2,sst[0:50]] medianlowstrp = median(lowstrp) help, medianlowstrp,fp_wl_upl if (first) then oldmed=medianlowstrp IF medianlowstrp LT 0.4 or medianlowstrp GT 2*oldmed THEN BEGIN ;IF 0 THEN BEGIN print,string(7b) ; bell toskip='y' print,'Encountered bad lowstrp.' ;wait,2 ;debugtime=5 ;IF toskip EQ 'y' THEN BEGIN IF 0 THEN BEGIN print,'Skipping ',filen wait,2 GOTO, skipped ENDIF ENDIF teststrip=reform(strp[ss,*]) stripends(0:10)=teststrip[0:10] stripends(11:21)=teststrip[ht-11:ht-1] newmedend = median(stripends) IF median(stripends(0:10)) LE 0 THEN newmedend = median(stripends(11:21)) IF median(stripends(11:21)) LE 0 THEN newmedend = median(stripends(0:10)) ; added 981116, NBR newmed = median(strp(ss,*)) IF ~first THEN oldmed = median(stripn0(ss,*)) IF first THEN printf,l1,obs_timen+' '+filename+' box_ref = '+trim(box_ref)+', rconv= '+trim(rconv) logmsg= obs_timen+' '+ filename+ $ ': carrmap_x ='+ STRING(carrmap_x,'(i3)')+ $ ', factor = ' + string(factor,'(f5.2)') + $ ', med(strip) = '+ string(newmed,'(f5.2)') + $ ', med(im) = '+ STRING(m,'(f5.2)') + $ ', max(im) = ' + STRING(mx,'(f5.2)')+ $ ', LP='+ string(lpnum,'(i2)')+ offinfo printf,l1, logmsg IF stdev(stripends) LT 0.0001 THEN BEGIN print,string(7b) ; bell print,'Stdev indicates Missing data. Skipping image?' wait,2 printf,l1,obs_timen+' '+filename+': Stdev indicates Missing data.' toskip='y' ENDIF IF newmed LE 0 THEN BEGIN print,string(7b) ; bell print,'Strip is mostly negative or zero; skipping image?' debugtime=30 ;wait,2 printf,l1,obs_timen+' '+filename+': Strip is mostly negative or zero.' toskip='y' ENDIF ; Test for strange blocks in image IF avg(abs(teststrip-median(teststrip))) GT 0.18 THEN BEGIN ; do not skip images where discontinuity is a few pixels print,string(7b) ; bell print,'Continuity test fails. Skipping image?' ;debugtime=30 ;wait,2 printf,l1,obs_timen+' '+filename+': Continuity test fails.' toskip='y' ENDIF w=where(oldmedends) ww =where(oldstripmeds) IF w(0) LT 0 THEN est_val = newmedend ELSE BEGIN ;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,'ends =',newmedend,', strip median =',newmed print,'End est_val =',est_val, ', Strip _est_val =',str_est_val strip_per = (newmed - str_est_val)/str_est_val ; Percent variance of strip from expected IF ( strip_per GT deltend) and not(keyword_set(RUNNING_DIFF)) and s(ord(k)).naxis2 GT 512 THEN toskip = 'y' ; changed from 3*str_est_val, 6/8/98 end_per = ABS(newmedend - est_val)/est_val ; Percent variance of end from expected IF cam EQ 'c1' THEN BEGIN IF newmed LT stripmin or newmed GT stripmax OR $ strip_per GT 0.8 or strip_per LT -0.3 THEN toskip = 'y' help,stripmin,stripmax,end_per,strip_per ENDIF ELSE BEGIN IF (( end_per GT (deltend+0.004))) or $ ; AND cn LT 1934 ( strip_per GT 0.4 or strip_per LT -0.15) THEN toskip = 'y' ; changed from 3*est_val, 6/8/98 IF median(oldmedends) EQ 0 and (medianok) THEN toskip='n' help,end_per,strip_per ENDELSE wait, nomwait ; for case where area of strip is masked and non-zero IF newmedend EQ avg(stripends) or newmed EQ m THEN toskip = 'y' ; nbr, 8/2/00 missing=where(strp le 0,nmissing) percentmissing=nmissing/(5.*181) IF nmissing GT 10 THEN BEGIN help,missing,percentmissing,debugtime message,'Replacing MISSING with '+trim(m),/info IF percentmissing GT maxmissing THEN BEGIN toskip='y' nooverride=1 ENDIF ELSE wait,1 ; 11% of strips is missing. strp[missing]=m ENDIF IF toskip NE 'y' and cx gt 1000 then stop IF toskip EQ 'y' and skipflag NE 'n' and skipwait LT 0 and (nooverride) THEN Begin print,string(7b) ;IF skiplast GT 5 AND mapdelta GT 5 THEN BEGIN ;toskip = 'n' ;print,string(7b) ;print,'Skiplast GT 5 AND mapdelta GT 5; not skipping' ;ENDIF ;IF (strip_per LT .8 or end_per LT .8) and not(keyword_set(RUNNING_DIFF)) or (carrmap_x LT 5) THEN BEGIN ;read,' Skip this one? (Default is yes) ',toskip ;IF toskip EQ '' THEN toskip = 'y' ;ENDIF ;wait,debugtime print,'To not skip, Ctrl-c and skipwait=1 ; =5 to not skip next 5 images.' print,'To continue, Ctrl-c and .cont' help,debugtime wait,debugtime IF (skipwait LT 1) THEN BEGIN ;print,string(7b) ; bell skiplast = skiplast+1 printf,l1,'Strip is too different; skipping image.' GOTO, skipped ENDIF ENDIF ;ELSE wait,2 skiplast = 0 ; ** END check strip ** 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 =',720-carrmap_x,' mapdelta =',mapdelta 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(1,*) = strp(r,*) tempstrip(0,*) = stripn0(r,*) IF mapdelta GT (mapsize/360)*15 or k EQ 0 or total(stripn0(r,*)) LE 0 THEN BEGIN ;current = cmap(*,*,r) ;no_zero = WHERE (current GT 0) ;IF no_zero(0) LT 0 THEN fill = 1.0 ELSE fill = median(current(no_zero)) tempstrip(*,*) = !values.f_nan ;fill ; ENDIF tempstrip(*,ht-barht:ht-1) = !values.f_nan ;max(cmap)>5.0 dummy = CONGRID(tempstrip,mapdelta*2,ht,/interp) newstrip = dummy(0:mapdelta,*) IF k GT 0 and NOT(first) 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,*) ;IF k EQ 0 and carrmap_x GT 0 THEN cmap(0,*,r)=strp(r,*) stripn0(r,*) = strp(r,*) ;if mapdelta gt 5 then stop ; ; ** End carrwarp2. ** wset,0 ENDFOR tvscl,cmap(*,*,ss)mapmin ;,0,(ht+10)*(r) 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) ;IF s(ord(k)).naxis2 GT 512 THEN oldstripmeds(9) = newmed ELSE oldstripmeds(9)=1.01 oldstripmeds[9] = newmed skiplast = 0 ; 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) 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(*,*) ;stop ENDFOR ; r loop wset,0 tvscl,cmap(*,*,ss)mapmin,0,185*(r) ENDIF IF cam EQ 'c1' THEN $ printf,l1,'Carrmap_x='+string(carrmap_x)+'; on: '+filename+' @ '+obs_timen+' - off: '+off_hdr.filename+' @ '+off_hdr.date_obs+' '+off_hdr.time_obs $ ELSE $ printf,l1,'Carrmap_x='+string(carrmap_x)+'; '+filename+' @ '+obs_timen ctr=ctr+1 obs_timen0 = obs_timen first = 0 medianhist = shift(medianhist,-1) medianhist[9]=m median0=median(medianhist[where(medianhist gt 0)]) skipped: ENDFOR ;k loop end_radius = radius end_center(0) = cx end_center(1) = cy IF cam eq 'c1' THEN $ rollcor[1]=get_roll_or_xy(hdr,'roll',/deg) ELSE $ rollcor[1]=rotangle dhdr = scc_update_history(hdr,'$Id: carrmapmaker2.pro,v 1.15 2017/12/19 23:01:47 nathan Exp $') dhdr = scc_update_history(dhdr,samplemethod) print,'time0 = ',time0_str,' time1 = ',time1_str print print,'Found ',numim,' images for CR ',cn, ' ',wlimb print print,'Fill =',fill printf,l1,'Fill =',fill ;mnm=lasco_readfits('/net/corona/cplex3/monthly/c2min.fts',hm) ;diffcarr,cmap,time1_str,cam,rad,limb free_lun,l1 l1=99 ; to fix problems with common block later on free_lun,l2 ;stop cmap = reverse(cmap) print END ; carrmapmaker2