;+ ; :Description: ; procedure to average and demodulate comp data ; ; :Params: ; date_dir day of year to process, in YYYYMMDD format ; wave_type name of line for set of wavelengths ; ; :Keywords: ; synoptic synoptic processes subset of images in the morning ; ; :Author: Tomczyk ; ; :Modifier: Sitongia ;- ;+ ; :Description: ; Describe the procedure. ; ; :Params: ; date_dir ; wave_type ; ; :Author: sitongia ;- pro demod, date_dir, wave_type, synoptic = synoptic common comp_paths, bias_dir, flat_dir, mask_dir, inventory_basedir, binary_dir, $ defered_file, hot_file, $ ldm_basedir, usb_basedir, raw_synoptic_basedir, raw_usb_basedir, process_synoptic_basedir, process_usb_basedir, $ archive_dir, movie_dir, fullres_dir, log_dir common comp_constants, nx, ny, $ center1074, center1079, center1083, $ stokes, n_stokes, $ debug, logFileUnit ;; ADDED BY JAMES MICHAELIS ;; ADD COMMON BLOCKS FOR EACH LOGGER common logBlock_artifact, Counter_1, NameRegistry_1, IDRegistry_1, LUN_1 common logBlock_process, Counter_2, NameRegistry_2, IDRegistry_2, LUN_2 common logBlock_activity, Counter_3, NameRegistry_3, IDRegistry_3, LUN_3 common logBlock_observation, Counter_4, NameRegistry_4, IDRegistry_4, LUN_4 common logBlock_entity, Counter_5, NameRegistry_5, IDRegistry_5, LUN_5 common logBlock_dataset, Counter_6, NameRegistry_6, IDRegistry_6, LUN_6 common logBlock_qualityAssertion, Counter_7, NameRegistry_7, IDRegistry_7, LUN_7 common logBlock_qualityEvidence, Counter_8, NameRegistry_8, IDRegistry_8, LUN_8 common logBlock_fitsHeader, Counter_9, LUN_9 ;; END ;; ADDED BY JAMES MICHAELIS ;; RESTORE VARIABLE STATE FROM PRIOR EXECUTION (allows logBlock values to be carried over) Restore,'CoMP_commonBlocks.sav' ;; END COMPILE_OPT IDL2 ; Configure comp_paths comp_initialize, date_dir printf, logFileUnit, 'demod' raw_dir = raw_usb_basedir + date_dir process_dir = process_usb_basedir + date_dir device,decomposed=0 ans=' ' cal='no' ;are these calibration images? (yes or no) stray='no' ;make stray light correction? (yes or no) avg_method='mean' ;averaging method (mean, median, min, medavg, minavg) skip_first = 'no' ; skip first image at each wavelength ('yes') or just first image in recipe ('no') use_flat_center = 'yes' ; Use the occulter center information from the flat file instead of finding it here beam_multiplies_wave = 'yes' ; multiply beam times wave to get unique waves infiles=process_dir+'/'+wave_type+'_files.txt' ;file with list of file names flatfile='flat.fts' ;file to use for flats extens=wave_type ;entension for files ; synoptic processing works on a subset of files in the morning n_file_lines = file_lines(infiles) if keyword_set(synoptic) then n_file_lines <= 40 ; compute transformation arrays for distortion removal x=rebin(findgen(nx),nx,nx) y=transpose(x) k1=0.99353 x1new=x*.5*(1.+k1)+y*.5*(1.-k1) y1new=x*.5*(1.-k1)+y*.5*(1.+k1) k2=1.00973 x2new=x*.5*(1.+k2)+y*.5*(1.-k2) y2new=x*.5*(1.-k2)+y*.5*(1.+k2) ; compute position of outside pixels for stray light removal if stray eq 'yes' then begin x1024=rebin(indgen(1024)-512.,1024,1024) y1024=transpose(x1024) outside=where(abs(y1024-x1024) gt 890 or abs(x1024+y1024) gt 670) endif ; open files openr,1,infiles ;file containing filenames openw,2,process_dir+'/demod.log' ;log file ; loop over filenames name=' ' for file_count = 0, n_file_lines - 1 do begin ; loop for each file ------------------------ readf,1,name,format='(a19)' printf, logFileUnit,name printf,2,name,' ',systime(/UTC) fits_open,raw_dir+"/"+name,fcb ;open input fits file num_images=fcb.nextend if debug eq 1 then print,num_images,' images in file' ; take inventory of images in this file inventory,fcb,beam,group,wave,pol,type,expose,cover,cal_pol,cal_ret ; identify unique wavelengths, polarizations and beams if debug eq 1 then print,'identify unique wavelengths, polarizations and beams' uniq_wave=wave[uniq(wave,sort(wave))] num_wave=n_elements(uniq_wave) uniq_pol=pol[uniq(pol,sort(pol))] num_pol=n_elements(uniq_pol) uniq_beam=beam[uniq(beam,sort(beam))] num_beam=n_elements(uniq_beam) if debug eq 1 then begin print,uniq_pol print,uniq_wave print,uniq_beam endif ; Do not use first image at each wavelength (put them in group -1) if debug eq 1 then print,'do not use first image at each wavelength (put them in group -1)' wave[0]=-1 ;do not use first image group[0]=-1 ;or do not use first image if skip_first eq 'yes' then begin bw = wave*float(beam) ;multiply wavelength by beam sign to allow to find uique wavelengths uniq_bw=bw[uniq(bw,sort(bw))] ;find unique wavelengths/beams if debug eq 1 then print,uniq_bw num_bw=n_elements(uniq_bw) for i=0,num_bw-1 do begin good=where(bw eq uniq_bw[i]) group[good[0]]=-1 endfor endif which0=where(group eq 0,num0) if num0 eq 0 then group=group-1 ;if no images in group 0 then decrement group if debug eq 1 then for i=0,n_elements(beam)-1 do $ print,beam[i],group[i],wave[i],pol[i],type,expose,cover,cal_pol ; read primary header and get time and exposure time if debug eq 1 then print,'read header and get time and exposure time' fits_read,fcb,d,primary_header,/header_only,exten_no=0 date_str=sxpar(primary_header,'DATE_OBS') time_str=sxpar(primary_header,'TIME_OBS') fits_read,fcb,d,header,/header_only,exten_no=1 exposure=sxpar(header,'EXPOSURE') ; interpret date and time if debug eq 1 then print,'interpret date and time' month=fix(strmid(name,4,2)) day=fix(strmid(name,6,2)) year=fix(strmid(name,0,4)) if year lt 2000 then year=year+2000 if strlen(time_str) eq 10 then begin hours=fix(strmid(time_str,0,1)) mins=fix(strmid(time_str,2,2)) secs=fix(strmid(time_str,5,2)) endif else begin hours=fix(strmid(time_str,0,2)) mins=fix(strmid(time_str,3,2)) secs=fix(strmid(time_str,6,2)) endelse if hours lt 7 then hours=hours+12 time=float(hours)+float(mins)/60.+float(secs)/3600. if debug eq 1 then print,time_str,hours,mins,secs ; compute solar ephemeris quantities from date and time jd=tojd(day,month,year,hours,mins,secs) ephem2,jd,sol_ra,sol_dec,b0,p_angle,semi_diam,sid_time,dist,xsun,ysun,zsun ; interpolate dark image for time of observation if debug eq 1 then print,'interpolate dark image for time of observation' dark=dark_interp(date_dir,time,exposure) if debug eq 1 then print,'read dark image at:',time,' with exposure of:',exposure ; read flat images for appropriate for these wavelengths and time if debug eq 1 then print,'read flat images for appropriate for these wavelengths and time' read_flats,date_dir,wave,beam,time,flat,flat_header,flat_waves,flat_names,flat_expose,file=flatfile flat=flat*expose/flat_expose ;modify for exposure times ; Do once when using centers from flats if use_flat_center eq 'yes' then begin xcent1 = sxpar(flat_header,'OXCNTER1') - nx/2 ycent1 = sxpar(flat_header,'OYCNTER1') - 1024 + nx/2 radius1 = sxpar(flat_header,'ORADIUS1') xcent2 = sxpar(flat_header,'OXCNTER2') - 1024 + nx/2 ycent2 = sxpar(flat_header,'OYCNTER2') - nx/2 radius2 = sxpar(flat_header,'ORADIUS2') c1 = [xcent1,ycent1,radius1] c2 = [xcent2,ycent2,radius2] endif ; Field position fxcent1 = sxpar(flat_header,'FXCNTER1') - nx/2 fycent1 = sxpar(flat_header,'FYCNTER1') - 1024 + nx/2 fradius1 = sxpar(flat_header,'FRADIUS1') fxcent2 = sxpar(flat_header,'FXCNTER2') - 1024 + nx/2 fycent2 = sxpar(flat_header,'FYCNTER2') - nx/2 fradius2 = sxpar(flat_header,'FRADIUS2') fc1 = [fxcent1,fycent1,fradius1] fc2 = [fxcent2,fycent2,fradius2] ; set up matrix for image rotation x0=float(nx)/2. y0=float(nx)/2. x=rebin( findgen(nx)-x0,nx,nx ) y=transpose(rebin( findgen(nx)-y0,nx,nx )) angle=p_angle+180. xp=x*cos(angle*!pi/180.) -y*sin(angle*!pi/180.) yp=x*sin(angle*!pi/180.) +y*cos(angle*!pi/180.) ; make averages of unique combinations as designated by group if debug eq 1 then print,'make averages of unique combinations as designated by group' printf,2,'Start Making Averages ',systime(/UTC) ngroup=max(group)+1 if debug eq 1 then print,ngroup,' groups' dat1=fltarr(nx,nx,ngroup,/nozero) dat2=fltarr(nx,nx,ngroup,/nozero) pos1=fltarr(3,ngroup) pos2=fltarr(3,ngroup) group_pol=strarr(ngroup) ;arrays to hold pol, wave and beam for each group group_wave=fltarr(ngroup) group_beam=intarr(ngroup) num_in_group=intarr(ngroup) ; Number of images averaged in each group for ig=0,ngroup-1 do begin good=where(group eq ig,num) ;locate images in this group data=fltarr(1024,1024,num) ;set up array for all images in this group group_pol[ig]=pol[good[0]] group_wave[ig]=wave[good[0]] group_beam[ig]=beam[good[0]] num_in_group[ig]=num printf, logFileUnit,format='(i3," images with",f10.3,2x,a," Beam:",i3,2x,20i4)', $ num,group_wave[ig],group_pol[ig],group_beam[ig],good ; read in all images in this group, dark correct and average if debug eq 1 then print,'read in all images in this group, dark correct and average' for i=0,num-1 do begin fits_read,fcb,dat,header,exten_no=good[i]+1 if sxpar(header,'DEMULT') eq 0 then dat=demultiplex(dat) dat=float(dat)-dark dat=fixrock(dat,0.030) dat=fix_image(dat) data[*,*,i]=dat endfor s = size(data) if s[0] eq 3 then begin avg=form_average(data,method=avg_method) endif else begin avg=data endelse ; optionally make stray light correction if debug eq 1 then print,'optionally make stray light correction' if stray eq 'yes' then begin degree=0 if degree eq 0 then avg=avg-median(avg[outside]) else begin xfit=x1024[outside] yfit=y1024[outside] zfit=avg[outside] coefs=surf_fit_xy(zfit,xfit,yfit,degree) avg=avg-eval_surf(coefs,x1024,y1024) endelse endif if debug eq 2 then begin ;display raw data ba=bytscl(avg,0,2000) tvwin,ba profiles,ba endif ; make flat correction if debug eq 1 then print,'make flat correction' if beam_multiplies_wave eq 'yes' then begin iflat=where( flat_waves eq group_wave[ig]*group_beam[ig], count) endif else begin iflat=where( flat_waves eq group_wave[ig], count) endelse if count le 0 then begin message, 'Failed to find flat to use' endif printf, logFileUnit,'flat used:',iflat[0], ' named ', flat_names[iflat[0]] avg=avg/flat[*,*,iflat[0]] avg=fix_hot(avg, hot_file) d1=extract1(avg) ;extract sub-arrays d2=extract2(avg) ; remove distortion d1=interpolate(d1,x1new,y1new,cubic=-0.5,missing=0.) d2=interpolate(d2,x2new,y2new,cubic=-0.5,missing=0.) ; find position of images if use_flat_center eq 'no' then begin c1=find_image(d1) c2=find_image(d2) endif pos1[*,ig]=c1 pos2[*,ig]=c2 dat1[*,*,ig]=d1 dat2[*,*,ig]=d2 wait,0.1 endfor fits_close,fcb ; compute mean position of images and standard deviations if debug eq 1 then print,'compute mean position of images and standard deviations' xcent1=median(pos1[0,*]) & xcent2=median(pos2[0,*]) ycent1=median(pos1[1,*]) & ycent2=median(pos2[1,*]) radius1=median(pos1[2,*]) & radius2=median(pos2[2,*]) xsig1=stdev(pos1[0,*]) & xsig2=stdev(pos2[0,*]) ysig1=stdev(pos1[1,*]) & ysig2=stdev(pos2[1,*]) rsig1=stdev(pos1[2,*]) & rsig2=stdev(pos2[2,*]) ; xsig1=0.0 & xsig2=0.0 ; ysig1=0.0 & ysig2=0.0 ; rsig1=0.0 & rsig2=0.0 printf, logFileUnit,xcent1,ycent1,radius1,xcent2,ycent2,radius2 printf,2,xcent1,ycent1,radius1,xcent2,ycent2,radius2 printf, logFileUnit,xsig1,ysig1,rsig1,xsig2,ysig2,rsig2 printf,2,xsig1,ysig1,rsig1,xsig2,ysig2,rsig2 ; compute image offsets xpp1=xp+x0+xcent1 ypp1=yp+y0+ycent1 xpp2=xp+x0+xcent2 ypp2=yp+y0+ycent2 ; Masking ; Occulter mask dmask1=disk_mask(c1[2], dx=c1[0], dy=c1[1]) dmask2=disk_mask(c2[2], dx=c2[0], dy=c2[1]) ; Field mask field_mask_1 = field_mask(fradius1, dx=fxcent1, dy=fycent1) field_mask_2 = field_mask(fradius2, dx=fxcent2, dy=fycent2) ; New field masks, slightly larger, to create larger overlap to mask overlap_mask_1 = field_mask(fradius1 + 10., dx=fxcent1, dy=fycent1) overlap_mask_2 = field_mask(fradius2 + 10., dx=fxcent2, dy=fycent2) ; Identify the overlap of images, from the field positions tmp_img=fltarr(1024,1024) tmp_img[0:nx-1,1024-nx:1024-1]=tmp_img[0:nx-1,1024-nx:1024-1] + overlap_mask_1 tmp_img[1024-nx:1024-1,0:nx-1]=tmp_img[1024-nx:1024-1,0:nx-1] + overlap_mask_2 overlap = where(tmp_img eq 2) tmp_img[overlap] = 0. ; Extract sub-images for masking overlap_mask_1 = extract1(tmp_img) ; note that these have not been shifted like the data will later be overlap_mask_2 = extract2(tmp_img) ; Combine masks mask1 = dmask1 * field_mask_1 * overlap_mask_1 mask2 = dmask2 * field_mask_2 * overlap_mask_2 ; Translate and rotate masks mask1 = interpolate(mask1,xpp1,ypp1,missing=0.) mask2 = interpolate(mask2,xpp2,ypp2,missing=0.) ; Post mask pmask = post_mask(p_angle - 2.5) ; Combine all masks into one to be superset for all images mask = mask1 * mask2 * pmask ; Translate, rotate images and reorganize from beam swapping signal=fltarr(nx,nx,ngroup,/nozero) ;array for signal for each group back=fltarr(nx,nx,ngroup,/nozero) ;array for background for each group for ig=0,ngroup-1 do begin d1=dat1[*,*,ig] d2=dat2[*,*,ig] d1=interpolate(d1,xpp1,ypp1,missing=0.,cubic=-0.5) ;translate and rotate images d2=interpolate(d2,xpp2,ypp2,missing=0.,cubic=-0.5) ; subtract background, unless these are cal images in which case add them if debug eq 1 then print,'extract corona and continuum' ; associate coronal and continuum images with signal and back arrays ; add signal and back if these are cal images if cal eq 'yes' then signal[*,*,ig]=d2+d1 else begin if group_beam[ig] eq 1 then begin signal[*,*,ig]=d2 back[*,*,ig]=d1 endif else begin signal[*,*,ig]=d1 back[*,*,ig]=d2 endelse endelse endfor s = size(back) avg_back = total(back, 3) / s[3] mean_back=median(avg_back[where(mask eq 1.)]) ;compute mean background level ; demodulate into Stokes parameters for unique wavelengths and beams printf,2,'Start Demodulation ',systime(/UTC) if debug eq 1 then print,'Start Demodulation ',systime(/UTC) stokes=['I','Q','U','V'] ; array keeping track of which parameters were observed at which wavelengths stokes_there=fltarr(4,num_wave) stokes_there[0,*]=1 ;assume intensity is there ; array for demodulated stokes parameters ; ordered by size, stokes, wavelength, beam demod=fltarr(nx,nx,4,num_wave,num_beam) ; TODO store separate backgrounds background=fltarr(nx,nx,num_wave) ; Array to count number of images going into average num_averaged = intarr(4,num_wave) num_back = intarr(num_wave) ; array for demodulated stokes parameters after beam differencing demod_diff=fltarr(nx,nx,4,num_wave) for ig=0,ngroup-1 do begin if debug eq 1 then print,group_pol[ig],group_wave[ig],group_beam[ig] stokes_param=strmid(group_pol[ig],2,1) ;get stokes parameter and sign stokes_sign=strmid(group_pol[ig],1,1) ;for this group which_wave= where( group_wave[ig] eq uniq_wave,count ) if count le 0 then message, 'Failed to find which_wave' which_beam= where( group_beam[ig] eq uniq_beam,count ) if count le 0 then message, 'Failed to find which_beam' which_stokes=where( stokes_param eq stokes, count ) if count le 0 then message, 'Failed to find which_stokes' stokes_there[which_stokes,which_wave]=1 if debug eq 1 then print,format='(4(i3,2x))',$ which_stokes,which_wave,which_beam ; 1083: don't subtract back if wave_type ne '1083' then begin sig_back=signal[*,*,ig]-back[*,*,ig] endif else begin sig_back=signal[*,*,ig] endelse ; add intensity demod[*,*,0,which_wave[0],which_beam[0]]= $ demod[*,*,0,which_wave[0],which_beam[0]]+sig_back/float(num_pol) num_averaged[0,which_wave[0]] += num_in_group[ig] ; add stokes parameter with proper sign if stokes_sign eq '+' then begin demod[*,*,which_stokes[0],which_wave[0],which_beam[0]]= $ demod[*,*,which_stokes[0],which_wave[0],which_beam[0]]+sig_back/2. endif else begin demod[*,*,which_stokes[0],which_wave[0],which_beam[0]]= $ demod[*,*,which_stokes[0],which_wave[0],which_beam[0]]-sig_back/2. endelse num_averaged[which_stokes[0],which_wave[0]] += num_in_group[ig] ; Background by wavelength background[*,*,which_wave[0]] += back[*,*,ig]/float(num_pol) num_back[which_wave[0]] += num_in_group[ig] endfor ; compute demodulated difference between beams temporally if num_beam eq 2 then begin for iw=0,num_wave-1 do begin for ip=0,3 do begin demod_diff[*,*,ip,iw]=(demod[*,*,ip,iw,1]+demod[*,*,ip,iw,0])/2. endfor background[*,*,iw] /= 2. endfor endif else begin for iw=0,num_wave-1 do for ip=0,3 do begin demod_diff[*,*,ip,iw]=demod[*,*,ip,iw,0] endfor endelse ; write Stokes parameters to output file printf,2,'Start Writing Images ',systime(/UTC) if debug eq 1 then print,'Start Writing Images ',systime(/UTC) fits_open,process_dir+"/"+ut_filename(strmid(name,0,15))+'.comp.'+extens+'.fts',fcbout,/write ; Get rid of all the blank comments sxdelpar,primary_header,'COMMENT' ; Save date of processing month_name = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', $ 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] today = SYSTIME(/UTC) month = strmid(today,4,3) month = where(month EQ month_name) + 1 month = month[0] month = string(month) month = strmid(month,10,2) day = strmid(today,8,2) year = strmid(today,20,4) datecal = year + '-' + month + '-' + day ; Change spaces to zeros for n=0,strlen(datecal) do begin char = strmid(datecal,n,1) if (char EQ ' ') then strput,datecal,'0',n endfor sxaddpar, primary_header,'OBJECT','corona', ' Coronal Emission', AFTER='LOCATION' ; Change the processing level sxaddpar, primary_header,'LEVEL','L1', ' Processing Level' sxaddpar, primary_header,'DATE-CAL',datecal, " Date of calibration processing" sxaddpar, primary_header, 'DEMODVER',"$Revision$", ' demod Software Version' sxaddpar, primary_header, 'NTUNES', num_wave, " Number of wavelength tunings", BEFORE='TNELNGTH' sxaddpar, primary_header,'BUNIT','MILLIONTHS', ' Millions of brightness of solar disk' sxaddpar, primary_header,'METHOD',avg_method, ' Averaging method' sxaddpar, primary_header, "COORDNAM", " HELIOCENTRIC", "COORDINATE SYSTEM NAME" sxaddpar, primary_header, "CTYPE1", "X", " AXIS 1 TYPE: X [EAST->WEST ] GEOCENTRIC" sxaddpar, primary_header, "CTYPE2", "Y", " AXIS 2 TYPE: Y [SOUTH->NORTH] GEOCENTRIC" sxaddpar, primary_header, "CRPIX1", nx/2+0.5, " X [EAST->WEST ] SUN CENTER [INDEX:PIXELS]" sxaddpar, primary_header, "CRPIX2", ny/2+0.5, " Y [SOUTH->NORTH] SUN CENTER [INDEX:PIXELS]" sxaddpar, primary_header, "CROTA1", 0.0, " X [EAST->WEST ] ROTATION [DEG.] FROM REFERENCE" sxaddpar, primary_header, "CROTA2", 0.0, " Y [SOUTH->NORTH] ROTATION [DEG.] FROM REFERENCE" sxaddpar, primary_header, "OCRAD1", radius1, " [pixels] median Occulter Radius sub-image 1" sxaddpar, primary_header, "OCRAD2", radius2, " [pixels] median Occulter Radius sub-image 2" sxaddpar, primary_header, "FCRAD1", fradius1, " [pixels] Field Radius sub-image 1" ; Center of field, offset from center of occulter and rotated by angle fxcent1p = nx / 2.0 + fxcent1*cos(angle*!pi/180.) - fycent1*sin(angle*!pi/180.) fycent1p = ny / 2.0 + fxcent1*sin(angle*!pi/180.) + fycent1*cos(angle*!pi/180.) sxaddpar, primary_header, "FCENX1", float(fxcent1p), " [pixels] Field center X sub-image 1" sxaddpar, primary_header, "FCENY1", float(fycent1p), " [pixels] Field center Y sub-image 1" sxaddpar, primary_header, "FCRAD2", fradius2, " [pixels] Field Radius sub-image 2" ; Center of field, offset from center of occulter and rotated by angle fxcent2p = nx / 2.0 + fxcent2*cos(angle*!pi/180.) - fycent2*sin(angle*!pi/180.) fycent2p = ny / 2.0 + fxcent2*sin(angle*!pi/180.) + fycent2*cos(angle*!pi/180.) sxaddpar, primary_header, "FCENX2", float(fxcent2p), " [pixels] Field center X sub-image 2" sxaddpar, primary_header, "FCENY2", float(fycent2p), " [pixels] Field center Y sub-image 2" sxaddpar, primary_header, "CRRADIUS", float(semi_diam)/4.46, " [pixels] Solar Radius" sxaddpar, primary_header, "CDELT1", 4.46, " solar_X coord increment [arcsec/pixel]" sxaddpar, primary_header, "CDELT2", 4.46, " solar_Y coord increment [arcsec/pixel]" ; Occulter ID and size occulter_id=sxpar(header,'OCCULTER') sxdelpar,primary_header,'OCCULTER' sxaddpar, primary_header, "OCC-ID", occulter_id, " Occulter Identification Number" occulter_size = occulter(occulter_id) sxaddpar, primary_header, "OCC-SIZE", occulter_size, " [mm] Occulter size" sxaddpar,primary_header,'BACKGRND',mean_back, " Mean background" sxaddpar,primary_header,'RSUN',float(semi_diam), " [arcsec] Solar Radius" sxaddpar,primary_header,'SOLAR_P0',float(p_angle), " [degrees] P Angle" sxaddpar,primary_header,'SOLAR_B0',float(b0), " [degrees] B Angle" ; TODO sxaddpar,header,'SOLAR_L0',float(l0), " [degrees] L Angle" sxaddpar,primary_header,'SOLAR_RA',float(sol_ra), " [HOURS] Solar Right Ascension" sxaddpar,primary_header,'SOLARDEC',float(sol_dec), " [degrees] Solar Declination" ; Fix the date/time in UT fix_header_time, primary_header sxdelpar,header,'COMMENT' sxdelpar,header,'SEQUENCE' sxdelpar,header,'BEAM' sxdelpar,header,'LCVR1VOL' sxdelpar,header,'LCVR2VOL' sxdelpar,header,'LCVR3VOL' sxdelpar,header,'LCVR4VOL' sxdelpar,header,'LCVR5VOL' sxdelpar,header,'LCVR6VOL' sxdelpar,header,'LCVR1TMP' sxdelpar,header,'LCVR2TMP' sxdelpar,header,'LCVR3TMP' sxdelpar,header,'LCVR4TMP' sxdelpar,header,'LCVR5TMP' sxdelpar,header,'LCVR6TMP' ; Write the input primary header into the output fits_write, fcbout, 0, primary_header for ip=0,3 do for iw=0,num_wave-1 do begin if stokes_there[ip,iw] eq 1 then begin sxaddpar,header,'POLSTATE',stokes[ip], " POLARIZATION STATE" sxaddpar,header,'WAVELENG',uniq_wave[iw], " [NM] WAVELENGTH OF OBS" sxaddpar,header,'NAVERAGE',num_averaged[ip,iw], " Number of images averaged together", AFTER="EXPOSURE" sxaddpar,header,'FLATFILE',flat_names[iw], " Name of flat field file" sxaddpar,header,'DATAMIN',min(mask*demod_diff[*,*,ip,iw])," MINIMUM DATA VALUE" sxaddpar,header,'DATAMAX',max(mask*demod_diff[*,*,ip,iw])," MAXIMUM DATA VALUE" ename=stokes[ip]+', '+string(format='(f7.2)',uniq_wave[iw]) ; TODO mask here fits_write,fcbout,mask*demod_diff[*,*,ip,iw],header,extname=ename endif end sxdelpar,header,'POLSTATE' sxdelpar,header,'FLATFILE' for iw=0,num_wave-1 do begin sxaddpar,header,'WAVELENG',uniq_wave[iw], " [NM] WAVELENGTH OF OBS" sxaddpar,header,'NAVERAGE',num_back[iw], " Number of images averaged together", AFTER="EXPOSURE" sxaddpar,header,'DATAMIN',min(mask*background[*,*,iw])," MINIMUM DATA VALUE" sxaddpar,header,'DATAMAX',max(mask*background[*,*,iw])," MAXIMUM DATA VALUE" fits_write,fcbout,mask*background[*,*,iw],header,extname='B, '+string(format='(f7.2)',uniq_wave[iw]) ;write individual background image endfor fits_close,fcbout flush, logFileUnit flush,2 ;; ADDED BY JAMES MICHAELIS ;; ADD LOG STATEMENTS newName = ut_filename(strmid(name,0,15))+'.comp.'+extens+'.fts.gz' ID_Process_1 = log_process('Demod','procedure to average and demodulate comp data.', $ 'demod.pro','https://subversion.ucar.edu/HAO/MLSO/CoMP/Pipeline/demod.pro','3798','Interactive Data Language (IDL)','Steve Tomczyk',[loggedID('artifact','dark.fts'), loggedID('artifact','flat.fts')],[],[loggedID('dataset','Files for 6/29')]) ID_Artifact_1 = log_artifact(newName,'CoMP image entry processed by demod routine','FITS','',[ID_Process_1],newName) ID_Activity_1 = log_activity('Flat Correction','Account for detector & optics bias',[ID_Process_1]); ID_Activity_2 = log_activity('Stokes Demodulation','Extract intensity images from background and signal entries in data files',[ID_Process_1]); ID_Activity_3 = log_activity('Image Rotation','Rotate images to account for angle of solar inclination',[ID_Process_1]); ID_Activity_4 = log_activity('Dark Correction','Correct for dark current',[ID_Process_1]); ID_Activity_5 = log_activity('Distortion Correction','Correct for perceived irregularity in shape of sun',[ID_Process_1]); ;; END endfor ; Compress files printf, logFileUnit, "Compressing FITS files." spawn, 'gzip -f ' + process_dir + '/*.comp.' + wave_type + '.fts' close,1 close,2 printf, logFileUnit,'done' free_lun, logFileUnit stop_logger ; closes LUNs for each log Save, /ALL, filename='CoMP_commonBlocks.sav' ; save all variables (including common blocks) to file to be opened later. this preserves logging blocks for upcoming IDL scripts ;; END end ; Main program for testing demod, '20110808', '1074', /synoptic ;demod, '20110730', '1074' end