;+ ; ; Project : s3drs reconstruction software ; ; Name : p_scchdr ; ; Purpose : part of the s3drs GUI ; ; Explanation: Converts many standard FITS headers (esp. SECCHI & ; LASCO) to our internal geometry and structures for reconstruction ; ; Use : called from 's3drs.pro' and 'quickcme.pro' ; ; Inputs : hdr = a FITS header in structure format ; (optional) pxwrapper = our geometry structure ; (optional) iele = element ID in that pxwrapper structure ; ; Outputs : (optional) srjr structure of the values ; ; Keywords : ; ; Common : none ; ; Restrictions: none ; ; Side effects: updates 'pxwrapper' with the desired values, if given. ; ; Category : gui, FITS ; ; Prev. Hist. : None. ; ; Written : Sandy Antunes, NRC, March-April 2009 ; ;- ; NEEDS SOHO polarization translation to SECCHI type notations. ; If given the optional 'pxwrapper' structure plus an 'iele' index, ; it will put 'srjr' into that wrapper at position 'iele'. pro p_scchdr,hdr,srjr,N=N,pxwrapper=pxwrapper,iele=iele au2km=149598000.0 au2m=au2km*1000.0 if (n_elements(iele) eq 0) then iele=0 au2meters = 1.5E11; ; best guess here at FITS file processing level if (tag_exist(hdr,'bunit')) then begin bunit=hdr.bunit level='0.5' ; first guess if (bunit eq 'MSB') then level='1.0' end else begin bunit='MSB' level='sim 1.0' end sat_name='unknown' tele_name='unknown' if (tag_exist(hdr,'obsrvtry')) then begin sat_name=hdr.obsrvtry ; e.g. STEREO_A, STEREO_B if (sat_name eq '') then sat_name=hdr.instrume; hack for roto tomo end else if (tag_exist(hdr,'telescop')) then begin sat_name=hdr.telescop; e.g. SOHO end else if (tag_exist(hdr,'program')) then begin sat_name=hdr.program end if (tag_exist(hdr,'instrume')) then begin tele_name=hdr.instrume; e.g. LASCO, SECCHI end else if (tag_exist(hdr,'program')) then begin tele_name=hdr.program end detector=hdr.detector; e.g. COR1, COR2, EUVI, HI1, HI1, C2, C3, EIT ; GET COORDINATES ; we now use a non-rotating coordinate system like HAE or HCI. ; for coordinates, get HEE x,y,z, convert to HAE as an ; absolute frame, then calculate spherical rho, phi, d ; and cylindrical z0. ; ; The reason we keep HAE as our x/y/z is that it allows us to ; compare across models at different times, because it is an ; absolute reference frame ; if not HEE, check for HEEQ ; if not HEE x,y,z, convert to HE x,y,z first ; HEE = Heliocentric Earth Ecliptic ; conversion is just, e.g. ; CONVERT_STEREO_COORD, '2007-05-06T11:30', COORD, 'HCI', 'HEE' if (tag_exist(hdr,'simdate')) then begin obsdate=hdr.simdate end else begin get_fits_time,hdr,obsdate end units='m'; default is units in meters if (tele_name eq 'LASCO') then begin o1 = get_orbit(obsdate) ; lasco has GCI, GSE, GSM, HEC, CAR tempxyz = [o1.gse_x,o1.gse_y,o1.gse_z] CONVERT_STEREO_COORD, obsdate, tempxyz, 'GSE', 'HAE' tempxyz = tempxyz * 1000.0; convert LASCO GSE km to our meters end else if (tele_name eq 'SECCHI') then begin tempxyz = [hdr.heex_obs,hdr.heey_obs,hdr.heez_obs] CONVERT_STEREO_COORD, obsdate, tempxyz, 'HEE', 'HAE' end else begin ; unknown or sim data if (tag_exist(hdr,'haex_obs')) then begin tempxyz = [hdr.haex_obs,hdr.haey_obs,hdr.haez_obs] if (tag_exist(hdr,'hae_unit')) then units=hdr.hae_unit if (units eq 'au' or units eq 'AU') then begin ; convert to meters tempxyz=tempxyz*au2meters units='m' end end else if (tag_exist(hdr,'heex_obs')) then begin tempxyz = [hdr.heex_obs,hdr.heey_obs,hdr.heez_obs] CONVERT_STEREO_COORD, obsdate, tempxyz, 'HEE', 'HAE' end else if (tag_exist(hdr,'crln_obs')) then begin x = hdr.crln_obs ; angle of satellite from x-axis in degrees y = hdr.crlt_obs ; angle of satellite from y-axis in degrees z=hdr.dsun_obs * au2meters ; satellite distance in AU tempxyz = [x, y, z] CONVERT_STEREO_COORD, obsdate, tempxyz, 'Carrington', 'HAE' end else begin print,'unknown coordinates, exiting.' stop end end hae_coords = ptr_new(tempxyz); save, handy for later plotting xyz2pixon,tempxyz,rho,phi,LL,z0,ecliptic x=tempxyz[0] y=tempxyz[1] z=tempxyz[2] ; Also set whether image is aligned to solar north or ecliptic north ; Used in p_process_filelist to determine whether images need ; to be rotated to solar north or not. ; default for sims, SECCHI 0.5, SECCHI 1.0, most others align='ecliptic' if (tag_exist(hdr,'crota')) then crota=hdr.crota else $ if (tag_exist(hdr,'crota1')) then crota=hdr.crota1 else crota=0.0 rotateme=0 if (tele_name eq 'LASCO') then begin DEFSYSV, "!DRADEG", 180.0d/!DPI DEFSYSV, "!DDTOR", !DPI/180.0d ; angle=festival_ecl_solar_north(tempxyz[0],tempxyz[1],tempxyz[2]) ; Routine taken from festival_ecl_solar_north and embedded in here ; so we can use Pixon even if Festival is not installed. r = SQRT(x^2.0d + y^2.0d + z^2.0d) ecl_lat = ASIN(z/r) ecl_lon = ATAN(y, x) beta = !dtor*82.746294d lambda = !dtor*345.59337d sx = COS(beta)*SIN(lambda - ecl_lon) sy = COS(ecl_lat)*sin(beta) - SIN(ecl_lat)*COS(lambda - ecl_lon) angle = ATAN(sx, sy) * !radeg align = 'solar' + ' , ecl angle = '+strtrim(string(angle),2) + $ ' , crota = ' + strtrim(string(crota),2) rotateme=1 end else if (tele_name eq 'SECCHI' and level eq '1.0' and $ crota eq 0.0) then begin ; odd case of SECCHI Level1.0 where they set ROTATE_ON during secchi_prep align='solar' rotateme=1 end ; NOW FIND CENTERS ; LASCO ; next line is a hack for 0.5 files lacking some details if (tag_exist(hdr,'time_obs') eq 0 and tag_exist(hdr,'mid_time')) then $ hdr=addrep_struct(hdr,'time_obs',hdr.mid_time) if (tag_exist(hdr,'CRPIX1')) then begin ; extract properly from FITS header ; note that FITS goes 1 to N using image centers, while ; classic Pixon goes 0 to (N-1) using lower corner, therefore ; physical image center assuming N is even of N/2,N/2 ; translates in classic Pixon to (N/2-1),(N/2-1) k0=[hdr.CRPIX1-1.0, hdr.CRPIX2-1.0] ; optical center [2] if (sat_name eq 'SOHO' or sat_name eq 'LASCO') then begin sunc=get_sun_center(hdr) print,'sun center is ',sunc if (hdr.CROTA1 EQ 180) then begin center=512 ; (1024-1)/2. delx = center - sunc.xcen dely = center - sunc.ycen sunc.xcen = center + delx sunc.ycen = center + dely del = [center,center] - k0 k0 = center + del end print,'sun center is ',sunc sunk0=[sunc.xcen,sunc.ycen] end else if (strmid(sat_name,0,6) eq 'STEREO') then begin sunc=scc_sun_center(hdr) sunk0=[sunc.xcen,sunc.ycen] end else begin crval=[hdr.CRVAL1, hdr.CRVAL2] ; offset in arcsec crval=crval/hdr.cdelt1; convert to pixels sunk0=k0-crval end ; print,'debug: used crpix1' end else if (sat_name eq 'SOHO' or sat_name eq 'LASCO') then begin ; some work-arounds if there is trouble parsing the SOHO data k0=t_param(detector,'OCCENTER') sunk0str=t_param(detector,'CENTER') sunk0=sunk0str.center end else begin ; if no hdr info, assume it's a center-pointing sun-centered frame k0=[(N/2)-0.5,(N/2)-0.5]; guess at optical center [2] sunk0=[(N/2)-0.5,(N/2)-0.5]; guess at sun center [2] end ; NOW GET, ROTATE (IF NEEDED), AND REBIN DATA ;; print,sat_name,'debug: need better than not SOHO or LASCO' ; checked roll, we should not adjust for it (based on comparison ; of using this technique with Feb 01 2007 images, indicating that ; roll is already being accounted for in image processing somewhere. ; if (sat_name ne 'SOHO' and sat_name ne 'LASCO') then begin ; if (hdr.obsrvtry eq 'STEREO_B') then ab='B' else ab='A' ; angle=0.0 - get_stereo_roll(hdr.date_obs,ab,SYSTEM='HAE') ; rotateme=1 ; print,'Angle is ',angle ; end if (rotateme) then begin tangle=angle ; reposition optical center (since rotation was around sun center) deltak0=k0-sunk0 newk0=[0,0] newk0[0] = cos(tangle*!dtor)*deltak0[0] - sin(tangle*!dtor)*deltak0[1] newk0[1] = sin(tangle*!dtor)*deltak0[0] + cos(tangle*!dtor)*deltak0[1] k0=sunk0+newk0 end Ndata=hdr.naxis1 if (n_elements(N) eq 0) then N=Ndata scaling=float(Ndata)/float(N) if (N ne Ndata) then begin ; worked out here in long form just so it is clear how we figure this ; basectr=[(Ndata-1)/2.0,(Ndata-1)/2.0] ; newctr=[(N-1)/2.0,(N-1)/2.0] ; deltactr=k0-basectr ; newdeltactr=(deltactr/scaling) ; newk0=newctr+newdeltactr ; deltactr=sunk0-basectr ; newdeltactr=(deltactr/scaling) ; newsunk0=newctr+newdeltactr ; k0=newk0 ; sunk0=newsunk0 ; now using Brian's improved way k0 = (k0 - scaling + 1) / scaling sunk0 = (sunk0 - scaling + 1) / scaling end ;;; NOW set all the structure tags from the header ;? lla0=; sun-center in image, was [], default was [1,1,1]*0.5*(N-1) ;? d = length of voxel edge in AU? def 1.0 if (n_elements(pxwrapper) ne 0) then begin ; get lla0 from pxwrapper if possible lla0=pxwrapper.sr.lla0 L0 = pxwrapper.sr.L0 chl=pxwrapper.sr.chl end else begin ; default is 'center of cube' lla0=[1,1,1]*0.5*(Ncube-1) L0 = 214.943 ; aka 1 au in solar radii = constant, do not change chl=2 end ;given cdelt1 = arcsec of 1 pixel and rsun = radius of sun in arcsec, ; then 1 pixel subtends cdelt1/rsun in solar radii... so the total ; FOV is just cdelt1/rsun * N_pixels (in solar radii) and our ; L0tau is just this value/N or, cdelt1/rsun ; If rebinning, cdelt1' = cdelt1 * scaling, e.g. rebinning to 1/4th ; the pixels has 4x the arcsec per pixel ; Note this means L0tau is 'rsun per pixel' so rsun in pixels is just 1/L0tau ; ;orig L0tau = 60/N (12/N C2, 60/N C3);distance subtended by pixel at L0 in rsun ; first choose a default platescale, then choose a nicer one from header fovrsun=p_instrspec('fov',detector) L0tau=(2.0*fovrsun)/double(N) if (tag_exist(hdr,'rsun')) then begin ; e.g. SECCHI if (hdr.rsun eq 0) then begin ; not properly made, so use hard-coded version ; C1=4rsun FOV, C2=14rsunFOV end else begin ; more accurate and robust calculation ; L0tau = pixel in arcsec over rsun in arcsec ;rsun=hdr.rsun rsun=959.63 ; Sun size in arcsec at L (aka 1AU) is absolute L0tau = scaling*hdr.cdelt1/rsun end end else if (sat_name eq 'SOHO' or sat_name eq 'LASCO') then begin ;rsunp=get_solar_radius(hdr,/pixel); returns solar radius in pixels rsunp=959.63 ; Sun size in arcsec at L (aka 1AU) is absolute L0tau = scaling/rsunp end else if (tag_exist(hdr,'dsun_obs')) then begin ; SIM or direct calc ; tan(rsun) = rsun/distance, rsun=6.96x10^5 km, 1AU=150x10^6 km, ; 1 degree = 3600 arcsec, hdr.dsun_obs is already in meters, so ;rsun_arcsec = atan(6.96E8/hdr.dsun_obs) ; L0tau = pixel in arcsec over rsun in arcsec rsun_arcsec=959.63 ; Sun size in arcsec at L (aka 1AU) is absolute L0tau = scaling*hdr.cdelt1/rsun_arcsec end tau=0.5*L0tau/L0; e.g. 30/N/214.943 (or 0 for pinhole) ; this will probably break for LASCO e.g. polar='Clear' ptype=datatype(hdr.polar) ; PIXON needs polarizations of: ; 'b' brightness ; 'p' plane ; 'r' radial ; 't' tangential ; If not 'b', it also needs the angle as 'fan' ; if (tele_name eq 'SECCHI') then begin ; SECCHI polar=hdr.polar fan=0; starting default if (polar eq 2001) then begin polar='r' end else if (polar eq 2002) then begin polar='t' end else if (polar eq 1001 or polar eq 1002) then begin polar='b'; total brightness or polarization brightness end else if (polar eq 1003 or polar eq 1004) then begin ; percent polarized or polarization angle-- cannot use print,'Problem, polarization is ',polar,',exiting.' return end else begin ; otherwise, it is likely an individual polarized angle fan=polar ;rectify polarization angles to HAE ; note that the ground calibration curves indicate that ; the 'fan=0' of STEREO_B is 60 degrees off from the ; 'fan=0' of STEREO_A. So we apply the +-90 rotation ; of the CCD as well as the 60 degree relative offset ; to get A&B in the same alignment (but with no indication ; available as to where '0' is relative to the CCD). Thus ; this makes the 2 detectors co-aligned in terms of polarization ; even though we have no absolute '0' for that alignment. if (sat_name eq 'STEREO_B') then fan=fan+90+60 else $ if (sat_name eq 'STEREO_A') then fan=fan-90 if (fan lt 0) then fan=fan+360 if (fan gt 360) then fan=fan-360 polar='p' end end else if (datatype(hdr.polar) eq 'STR') then begin ; LASCO and some sims polar=hdr.polar if (polar eq 'Clear' or polar eq 'b' or polar eq 'B') then $ polar='b' else polar='p' ; need further work on LASCO polarizations here if (tag_exist(hdr,'fan')) then begin fan=hdr.fan ; e.g. 0, filter angle for plane-polarized frames end else begin fan=0 ; guess that angle is 0 if not otherwise knowable print,'Guessing polarization angle is 0, type is ',polar end end else begin ; generic or SIM fan=hdr.polar ; e.g. 0, filter angle for plane-polarized frames if (fan ne 0) then polar='p' else polar='b' print,'Guessing polarization angle is ',fan,', type is ',polar end ; polar= e.g. 'b', 'p', 't' or 'r' ; make masks for coronographs ; mask=congrid(quickC2mask(fname),N,N) ; (eps = pointing angle from sun center) imodel = 10; instrument-specific model, 10=B/B0 ; sets 'sr.an, sr.limbdark' as well as values needed for get_Istd_nd dvoxfrac=pxwrapper.sr.dvox/L0tau ; chl=-1; set chl to -1 for test/debug with large images o_flimb,chl,N,1.0/L0tau,dvoxfrac,lla0,fa,fb,fc,rmin2,an,limbdark Istd = get_istd_Nd(imodel,L0,1.0,L0tau,ImI0 = limbdark) ; date_obs now uses 'obsdate' from header, rather than the switch below ; need exptime in seconds, note level 1 tB uses '-1' but keeps ; the original exptime in expcmd if (tag_exist(hdr,'exptime')) then exptime=hdr.exptime else exptime=6.0 if (exptime eq -1) then if (tag_exist(hdr,'expcmd')) then $ exptime=hdr.expcmd else exptime=6.0 MSBperDN=get_msb2dn(hdr) if (sat_name eq 'SOHO' or sat_name eq 'LASCO') then DNperPhoton = 13.0 else $ DNperPhoton = 15.29 ; SECCHI value, also placeholder or best guess if (tag_exist(hdr,'biasdev')) then biasdev=hdr.biasdev else biasdev=0.0 ; set some placeholder/defaults for now eps=double(0) set = sin(eps)/tau roco = long(0); only for opixon as = long(1); ; z0=0.0; important, as z0 not properly implemented in cpixon srjr = { $ Ndata: N, $ scaling: scaling, $ ; tracks how original data was resampled sat_name: sat_name, $ tele_name: tele_name, $ detector: detector, $ L0tau: L0tau, $ tau: tau, $ rho: rho, $ ; same as gamma, but better descriptor gamma: rho, $ ; old pixon definitions, kept for legacy apps phi: phi, $ ; same as theta, but better descriptor theta: phi, $ ; old pixon definitions, kept for legacy apps LL: LL, $ z0: z0, $ polar: polar, $ fan: fan, $ k0: ptr_new(k0), $ sunk0: ptr_new(sunk0), $ Istd: Istd, $ limbdark: limbdark, $ imodel: imodel, $ date_obs: obsdate, $ MSBperDN: MSBperDN, $ DNperPhoton: DNperPhoton, $ biasdev: biasdev, $ hae_coords: hae_coords, $ exptime: exptime, $ ecliptic: ecliptic, $ align: align, $ roco: roco, $ eps: eps, $ set: set, $ as: as, $ level: level, $ bunit: bunit $ } ; optionally updates 'pxwrapper' if provided if (n_elements(pxwrapper) ne 0 and n_elements(iele) ne 0) then begin pxwrapper_update,pxwrapper,'sr',srjr,iele pxwrapper.sr.an=an; being a dbl(3), requires manual update end end