pro stksimages_sbsp,l1datadir, stks_data, outdir=outdir, no_save=no_save ; routine to generate quick-look images from Level1 SP data. A number ; of parameters are calculated from the spectral data that indicate ; apparent flux densities, polarization levels, field azimuth, etc. ; NOTE: this routine uses data an computations specific to the Hinode ; SP. The procedures have some hardwired parameters specific ; to that instrument that would need to be modified ; to be applicable to data from other instruments. ; INPUTS: l1datadir = path to directory containing the Level1 ; calibrated data ; outputs are FITS files, dimensioned [nx,ny] ([slit scan,along slit]) ; unless otherwise noted ; OUTPUTS: ; stks_data - structure containing 'savefile' equivilents, one tag per param ; no_save - remove the save file assumes return parameter 'stks_data' is supplied! ; avctr = empirical location of 6301.5 line center, pixels ; blapp = longitudinal apparent flux density, Mx cm^-2 ; btapp = transverse apparent flux density, Mx cm^-2 ; conti = continuum intensity, data units ; ctrline = 3D array [nx,ny,2]: pixel of minimum intensity for ; 6301.5, 6302.5 A lines ; dclntn = solar declination, arcseconds north of disk center ; dopcv = orbital Doppler shift DOP_CVR, pixels ; muval = cosine of heliocentric angle ; fitav = smoothed empirical spectral drift correction, 6301.5, = ; 29.0 - smoothed empirical spectral drift, pixels ; fitsp = smoothed predicted spectral drift correction, 6301.5, = ; 29.0 - smoothed predicted spectral drift, pixels ; fitww = smoothed empirical thermal drift along slit, pixels ; ftime = time from start of operation, minutes ; pii = Stokes I integrated over lines/conti ; plc = sqrt(Q^2+U^2) integrated over CONTINUUM/conti ; pll = sqrt(Q^2+U^2) integrated over lines/conti ; polnet = measure of linear polarization: (preferred-frame ; Stokes Q signal multiplied by qmask)/conti ; pq = Stokes Q integrated over lines/conti ; ptot = sqrt(Q^2+U^2+V^2) integrated over lines/conti ; pu = Stokes U integrated over lines/conti ; pv = Stokes V integrated over lines/conti ; rgtasn = solar right ascension, arcsec W of disk center ; rufaz = rough guess for field azimuth angle -- when ; linear pol. signals rotated to this angle, the ; signals are in the preferred frame. Radians ; CCL from solar west. ; scatlp = scattered light profile ; slat = solar latitude (carrington) degrees ; slong = solar longitude (carrington) degrees (from CM) ; sltdr = Ichimoto prediction, thermal drift along slit, pixels ; spcdr = Ichimoto/Kubo prediction, thermal spectral drift, pixels ; utimed = universal time from start of day, fractional hours ; vcs = 3D array [nx,ny,2]: pixel of Stokes V zero-crossing ; for 6301.5, 6302.5 A lines (has problems with ; anomalous profiles) ; wdelw = empirical thermal drift along slit, pixels ; History: ; 29-Aug-2007 - Bruce Lites ; 8-oct-2007 - S.L.Freeland - add OUTDIR & NO_SAVE keywords and stks_data output parameter ; 1-dec-2007 - Bruce Lites - compute position parameters (rgtasn, dclntn, ; muval, slat, slong, utimed) and scatt. light profile (scatlp) ; 23-jan-2008- Bruce Lites - simplify output save file by eliminating unneeded ; thermal drift variables (full set still in thermdsbsp.save) ; 25-jan-2008- Bruce Lites - use concat_dir to construct restore path names ; 14-jan-2009- Bruce Lites - if few or no low polarization pixels, construct ; the scat. light profile from lowest 64 pol. profiles ; include all the seldom-changed fixed parameters scal=get_logenv('SOT_SP_CALIBRATION') if scal eq '' then set_logenv,'SOT_SP_CALIBRATION',$ concat_dir('$SSW_SOT','calibration') no_save=keyword_set(no_save) if no_save and n_params() le 1 then begin box_message,["/NO_SAVE requested yet you did not supply an output parameter; I won't waste your time...", $ 'IDL> stksimage_sbsp,l1dir,stks_data,/no_save ; < return via stks_data'] return endif ; restore the calibration data polccc=file_search(get_logenv('$SOT_SP_CALIBRATION'),'polcal*geny') polccc=last_nelem(polccc) ; might want different logic in future? ; ; e.g. closest rather than latest if not file_exist(polccc) then begin box_message,'Cannot find polcal file under $SOT_SP_CALIBRATION' stop endif restgenx,file=polccc(0),tbbb,tpolnet,lbbb,lpv,qmask ; qmask is the Quiet-Sun Stokes Q profile mask. Note: this is not ; very good for large field strengths or inclined fields where MO ; effects can produce significant Stokes U and diminish the Q signal ; accordingly. BEWARE! ; get list of fits files in current directory derf = file_search(l1datadir , '*SP3D*.fits') ; tweak, S.L.Freeland 5-oct-2007 nfile = sizeof_sbsp(derf) print,'number of files = ',nfile ; extract the file names from the full path names fnames = strarr(nfile) ; establish output image array names for ifil = 0,nfile-1 do begin fnn = derf(ifil) spos = strpos(fnn,'SP4D20') fnames(ifil) = strmid(fnn,spos) endfor ; pixel ranges for integration of the Stokes signals ilo = intarr(4) & ihi = intarr(4) ; HARDWIRED PARAMETERS SPECIFIC TO HINODE SPECTRO-POLARIMETER ; calibration factor for pvpv as derived from inversions. Multiply ; pvpv by this number to get Mx cm^-2 ; NOTE: This is not currently used. Now we use the nonlinear ; calibration curve based upon ME data. This constant could be used ; for more accurate results when it is known that all fields are kG ; in strength, and that there are small fill fractions. If this ; constant is used on fields that are intrinsically weak, the results ; overestimate the apparent longitudinal flux by 30% cslp = 16836.5 ; setup parameters for PVPV computations nmult = 2 ;lctr = [27,73] ; Hardwired pixel centers for Hinode data. sp_prep.pro shifts ; the data to these pixel centers in producing L1 data. lctr = [29,75] ; Hardwired dispersion for Hinode SP, mA/pixel disper = 21.549 ; +/-300 mA wavelength range around spectral lines for integration of ; polarization signals liml = intarr(2,2) del = round(300./disper) for im = 0,nmult-1 do liml(*,im) = [lctr(im)-del,lctr(im)+del] ; continuum wavelength range ; avoid end of spectral range where spectrum is extended limc = [94,104] nwcont = limc(1) - limc(0) + 1 ; number of continuum wavelengths ilo(0) = limc(0) ihi(0) = limc(1) ; number of spectral pixels in spectral integration delnorm2 = nmult*(2*del+1) ; compute number of spectral pixels in spectral integration delnorm = 0 for im = 0,nmult-1 do begin delnorm = delnorm + liml(1,im) - liml(0,im) + 1 endfor ; limit for net polarization below which sum the scattered light profile pipmin = 0.0035 ; rotation angle of SP slit clockwise wrt Gband vertical dimension ; in degrees ;fgsprot = 0.2636 ; Now in header, don't hardwire fgsprot ; guesses for line center for azimuth estimation ceni = lctr(1) ; counter for number of averaged profiles for scattered light nscat = long(0) ks = -1 for kk = 0,nfile-1 do begin ks = ks+1 ; print,' start processing file ',ks fitsn = derf(kk) ; read and adjust level1 data for overflow in Stokes I, bitshifts readl1_sbsp,fitsn,dat,hdr ;dat = float(readfits(fitsn,hdr)) if ks eq 0 then begin ; store the first header of the sequence hdr0 = hdr nw = sizeof_sbsp(dat,1) & nsl = sizeof_sbsp(dat,2) ; temporary array for single stokes spectrum stks = fltarr(nw,4) ; output arrays stksimg = fltarr(nfile,nsl,4) rufaz = fltarr(nfile,nsl) polcont = stksimg rmsq = fltarr(nfile) rmsu = rmsq & rmsv = rmsq & rmsi = rmsq ; arrays for Stokes V zero-crossing, intensity minimum, ; continuum intensity, pvpv for each line vcs = fltarr(nfile,nsl,2) ctrline = fltarr(nfile,nsl,2) conti = fltarr(nfile,nsl) pv = fltarr(nfile,nsl) polnet = fltarr(nfile,nsl) pq = fltarr(nfile,nsl) pu = fltarr(nfile,nsl) pii = fltarr(nfile,nsl) ; wavelength average continuum net linear polarization signal plc = fltarr(nfile,nsl) ; wavelength average line net linear polarization sqrt(Q^2+U^2)/Ic pll = fltarr(nfile,nsl) ; wavelength average line total polarization sqrt(Q^2+U^2+V^2)/Ic ptot = fltarr(nfile,nsl) ; dclntn = solar declination, arcseconds north of disk center dclntn = fltarr(nfile,nsl) ; muval = cosine of heliocentric angle muval = fltarr(nfile,nsl) ; rgtasn = solar right ascension, arcsec W of disk center rgtasn = fltarr(nfile,nsl) ; scatlp = scattered light profile scatlp = dblarr(nw) ; slat = solar latitude (carrington) degrees slat = fltarr(nfile,nsl) ; slong = solar longitude (carrington) degrees (from CM) slong = fltarr(nfile,nsl) ; utimed = universal time from start of day, fractional hours utimed = fltarr(nfile,nsl) endif ; calculate the Stokes V zero crossing, Intensity minimum pixel, ; PVPV, continuum intensity for jj = 0,nsl-1 do begin ; if there are off-limb points identified by small intensities, don't ; try to fit. This will load zeros into the arrays and will produce ; some warning signals. Pick intensity of 500 as a safe low value, ; intensity in umbrae is larger than this. if mean(dat(1:10,jj,0) gt 500.,/double) then begin stks(*,*) = dat(*,jj,*) pvpvcomp_sbsp,stks,nmult,liml,limc,lctr,qmask,vcross, $ ctrlin,contin,pvpv,plpl,pqpq,pupu,pipi,plt,ptt,azmth ctrline(ks,jj,*) = ctrlin(*) vcs(ks,jj,*) = vcross(*) conti(ks,jj) = contin pv(ks,jj) = pvpv polnet(ks,jj) = plpl pu(ks,jj) = pupu pq(ks,jj) = pqpq pii(ks,jj) = pipi pll(ks,jj) = plt ptot(ks,jj) = ptt rufaz(ks,jj) = azmth endif endfor ; sum in the scattered light profile whrs = where(ptot(ks,*) le pipmin,countsct) if countsct gt 0 then begin for jj = 0,countsct-1 do begin indx = whrs(jj) scatlp = scatlp + dat(*,indx,0) endfor nscat = nscat + countsct endif ; sometimes the azimuth selection centers on -Q. Thus, for cases ; where polnet is negative, move the azimuth by 90 degrees and reverse ; the sign of polnet whrpp = where(polnet lt 0.,countp) pi2 = !pi/2. if countp gt 0 then begin tempruf = rufaz(whrpp) + pi2 whr90 = where(tempruf gt pi2,count90) if count90 gt 0 then tempruf(whr90) = tempruf(whr90) - !pi rufaz(whrpp) = tempruf polnet(whrpp) = -polnet(whrpp) endif ; do calculation of continuum net linear polarization ; veclin = sqrt(dat(*,*,1)^2 + dat(*,*,2)^2) ; average spectrally the sum of the squares, just as was ; done for the line polarization veclin = dat(*,*,1)^2 + dat(*,*,2)^2 for jj = 0,nsl-1 do plc(ks,jj) = $ sqrt(mean(veclin(ilo(0):ihi(0),jj),/double))/conti(ks,jj) ; ; UTIME COMPUTATION ; The output utimed is the universal time in hours. Here is how it ; is calculated in IDL from the header values of the Level1 FITS file: ; File clock time strtdate = ' ' strtdate = sxpar(hdr,'DATE_OBS') iyr = fix(strmid(strtdate,0,4)) imnth = fix(strmid(strtdate,5,2)) iday = fix(strmid(strtdate,8,2)) ahr = float(strmid(strtdate,11,2)) amn = float(strmid(strtdate,14,2)) asec = float(strmid(strtdate,17,6)) utimed(kk,*) = ahr + amn/60. + asec/3600. ; RGTASN COMPUTATION ; Right ascension, Arcseconds west of disk center. ; For Hinode, the spacecraft roll angle is such that the slit is ; nearly aligned N-S. But we must correct for the angle of the ; slit relative to the G-band image (fgsprot, to be SP Level1 ; keyword SPSLROT, but hardwired above), and also for the angle of the ; SP relative to solar north (CROTA2). Both angles measured ; clockwise from solar N. crota2 = sxpar(hdr,'CROTA2 ') fgsprot = sxpar(hdr,'FGSPROT ') ; sine of total rotation angle of slit sinrot = sin(!dtor*(crota2+fgsprot)) ; arcsec/pixel along the slit yscale = sxpar(hdr,'YSCALE ') ; solar right ascension of center of slit xcen = sxpar(hdr,'XCEN ') ; Range of pixels along the slit ; Note that the pixel indices along the slit direction are reversed from ; those during pre-launch testing. The images are now reversed in ; the vertical direction during reformatting. Also, pixel indices in ; the header like SPCCDIX0, SPCCDIX1 refer to the pre-launch CCD coordinates, ; not the coordinates of LEVEL0 data. spccdix0 = sxpar(hdr,'SPCCDIX0') spccdix1 = sxpar(hdr,'SPCCDIX1') ; define these pixel positions in terms of the LEVEL0 image coordinates ; reverse order so that low limit of active range is in spccd0 spccd0 = 1023-spccdix1 spccd1 = 1023-spccdix0 ; Extract the number of pixels summed along the slit (along serial direction) ssum = sxpar(hdr,'CAMSSUM') ; calculate the CCD pixel positions relative to the center of the ; CCD (pixel 511.5) ssll = float(spccd0 + indgen(nsl)*ssum) - 511.5 ; compute the right ascension for rotated slit rgtasn(kk,*) = xcen + sinrot*ssll*yscale/ssum ; DECLINATION COMPUTATION ; Declination, Arcseconds north of disk center. ; Center position of Y-range in Arcsec from disk center ycen = sxpar(hdr,'YCEN ') ; cosine of total rotation angle of slit cosrot = cos(!dtor*(crota2+fgsprot)) ; declination for digitized positions along slit dclntn(kk,*) = ycen + cosrot*ssll*yscale/ssum ; LATITUDE, LONGITUDE, MU COMPUTATION ; get the solar ephemeris values for this date/time. This gives ; the needed solar radius, P, and B0 angles. This routine also ; exists in FORTRAN77 as part of the ASP code, and was adapted ; from that code (a.aspcar.f) ; start loop over slit length for isl = 0,nsl-1 do begin sunframe_sbsp, iyr, imnth, iday, utimed(kk,isl), $ rgtasn(kk,isl), dclntn(kk,isl), $ sunlong, sunlat, valuemu, $ ras, dcs, gsdt, b0ang, peeang, ctrlong, r0 muval(kk,isl) = valuemu slong(kk,isl) = sunlong slat(kk,isl) = sunlat ; end loop over slit length endfor ; end loop over files endfor ; convert latitude, longitude to degrees slong = slong/!dtor slat = slat/!dtor ; if very few points below the threshold, find average profile of ; lowlmt profiles with lowest polarization lowlmt = 64 ; number of low-polarization profiles to average if nscat le lowlmt then begin print,' WARNING: number of points with pipmin<',pipmin,' < ',lowlmt print,' WARNING: using average profile of lowest pol. ',lowlmt,' pts' scatlp(*) = 0. psindx = sort(ptot) psindxl = psindx(0:lowlmt-1) ; find indices for these points ;mapindx = psindxl/nfile sltindx = psindxl/nfile ; find slit indices ;sltindx = psindxl - mapindx*nfile mapindx = psindxl - sltindx*nfile ; read through files again to average the low polarization profiles for kk = 0,lowlmt-1 do begin kmdx = mapindx(kk) ksdx = sltindx(kk) fitsn = derf(kmdx) ; read and adjust level1 data for overflow in Stokes I, bitshifts readl1_sbsp,fitsn,dat,hdr scatlp = scatlp + dat(*,ksdx,0) endfor nscat = lowlmt endif ; renormalize scattered light profile scatlp = scatlp/float(nscat) ; compute the apparent flux densities from the polarization signals bapp_sbsp,polnet,pv,tbbb,tpolnet,lbbb,lpv,btapp,blapp ; output the FITS files for each of the output parameters ; retain most of the Level1 header parameters ; image dimensions must be modified ; Should use hdr0 -- the first header in the sequence of files, ; as that contains the start times, etc. case 1 of data_chk(outdir,/string): ; user supplied no_save: outdir=get_temp_dir() else: outdir=curdir() endcase ; get the thermal drift parameters to add to the save file restore,concat_dir(l1datadir,'thermdrift.save') ; save results as a standard IDL-save file savename=concat_dir(outdir,'stksimg.save') save,filename=savename,blapp,btapp,conti,ctrline,dclntn, $ muval,pii,plc,pll,polnet,pq,ptot,pu,pv,rgtasn,rufaz,scatlp, $ slat,slong,utimed,vcs, $ ; thermal drift parameters appended wdelw,ftime,fitww,sltdr, $ spcdr,dopcv,fitsp,avctr,fitav if n_params() eq 2 then stks_data=ssw_save2struct(savename) if no_save then ssw_file_delete,savename return end