pro azam_hmi2afiles,outdir ;OUTPUT: outdir = path to output directory ; routine to generate asp-like a_* files for input to ; azam.pro from HMI inversion results ; common area for inversion results read in rdhmifits.pro common com_invertresults, $ fld, psi, azm, eta0, dmp, dop, bzero, b1mu, cen1, alpha, $ cct, chisq, occt, fld_er, psi_er, azm_er, $ cen1_er, alpha_er, xx, yy, slat, slong, $ yr, mnth, dy, hr, minit, sec, cdelt1, cdelt2,rconv ; get the array sizes sz = size(cct) nxx = sz(1) ; cct -> continuum intensity nyy = sz(2) ; for HMI, generate 1-D arrays of pixel rows for time ; variables hyr = intarr(nxx) hmnth = hyr & hdy = hyr & hhr = hyr & hminit = hyr hsec = fltarr(nxx) hyr(*) = yr & hmnth(*) = mnth & hdy(*) = dy hhr(*) = hr & hminit(*) = minit & hsec(*) = sec ;tim1 = systime(/seconds) ; ------------------- ; find local frame inclination (el0), azimuth (az0) print,' converting azimuth, inclination to local frame' print,' have patience, this may take a few minutes for large HMI images' azam_convertolocal,azm,psi,hyr,hmnth,hdy,hhr,hminit,hsec, $ slat,slong,rconv,az0,el0,az1,el1 ; ------------------- ;tim2 = systime(/seconds) ;print,'azam_convertolocal took ',tim2-tim1,' seconds' ; get solar ephemeris values for the start time for the header utim0 = hr(0) + minit(0)/60. + sec(0)/3600. solareph_sbsp,yr(0),mnth(0),dy(0),utim0,ras,dcs,gsdt,b0ang, $ peeang,ctrlong,r0 ; now read in the sample ASP header, modify the values, ; and write out again. Now hardwire the input header data ; so that we do not need an input file ;hdr = read_floats('sample_a___header') hdr = fltarr(128) hdr(*) = [6.0, 24.0, 92.0, 13.0, 256.0, 2.0, 11111.0, 11111.0, 1.0, $ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, $ 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, $ 0.0, 0.0, 0.0, 0.0, 0.0, 0.01, 1.0, 1.0, 0.1, $ 123.1, 45.3, 0.0, 12.59, 6302.50, 6301.51, 0.0, -300.0, 130.0, $ -300.0, 300.0, 0.0, -1.0, 4.0, 0.0, 0.0, 0.0, 0.0, $ 0.0, 0.0, -1.0, 300.0, 0.0, 0.0, 0.0, 0.0, 0.0, $ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, $ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, $ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, $ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 93.6201, 23.3985, 134.075, $ 2.19489, -5.47223, 257.435, 945.600, 182.977, -64.91, 0.0, 0.0, 0.0, $ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, $ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, $ 0.0, 0.0 ] dpr = 1./!dtor hdr(0) = float(mnth(0)) hdr(1) = float(dy(0)) hdr(2) = float(yr(0)) hdr(57) = float(nyy) hdr(96) = ras*dpr hdr(97) = dcs*dpr hdr(98) = gsdt*dpr hdr(99) = b0ang*dpr hdr(100) = peeang*dpr hdr(101) = ctrlong*dpr hdr(102) = r0 ; FOR NEW AZAM, determine the X-, Y-pixel spacing in arcseconds and ; write them into the AZAM a___header ; find values for arcseconds along x,y dimensions and insert them ; into the header for setting scale for plots. Use both xx,yy in ; case that observations are not oriented N-S, E-W ; range along x-dimension delx = sqrt( (xx(nxx-1,nyy/2)-xx(0,nyy/2))^2 +(yy(nxx-1,nyy/2)-yy(0,nyy/2))^2 ) delx = delx/float(nxx) dely = sqrt( (xx(nxx/2,nyy-1)-xx(nxx/2,0))^2 +(yy(nxx/2,nyy-1)-yy(nxx/2,0))^2 ) dely = dely/float(nyy) ; write them into apparently unused header locations 110,111 hdr(110) = delx ; arcseconds per pixel along X-direction hdr(111) = dely ; arcseconds per pixel along Y-direction print,'arcsec per pixel, X,Y: ',delx,dely ; write out the modified header err = write_floats(outdir+'a___header',hdr) ; now write out the a_* files for this op ; generate the a__points file temp = cct & temp(*,*) = 0. for ii = 0,nxx-1 do for jj = 0,nyy-1 do begin temp(ii,jj) = float(ii)*float(nyy) + float(jj) endfor err = write_floats(outdir+'a__points',temp) ; now set to zero the unsolved points whr = where(rconv eq 0,count) if count gt 0 then temp(whr) = 0. err = write_floats(outdir+'a_solved',temp) ; write the scattered light fraction in ASP form ; convert from fill fraction to the ASP scattered light fraction temp(*,*) = 1. - alpha err = write_floats(outdir+'a_alpha',temp) ; universal time ;for ii = 0,nxx-1 do temp(ii,*) = hr(ii) + minit(ii)/60. + sec(ii)/3600. for ii = 0,nxx-1 do temp(ii,*) = hr + minit/60. + sec/3600. err = write_floats(outdir+'a__utime',temp) ; write out doppler shifts, not line center positions, here err = write_floats(outdir+'a_cen1',cen1) ;err = write_floats(outdir+'a_cen2',cen2) ; net polarization ; !! NOTE !! old ASP considered PIP a percentage. Here it is ; a fractional value. This is taken into account in azam_op.pro ;err = write_floats(outdir+'a__pip',netpol) ; solar mu-value ;temp = cos(asin(sqrt(warc*warc+narc*narc)/float(r0))) temp = cos(asin(sqrt(xx*xx+yy*yy)/float(r0))) err = write_floats(outdir+'a__mu',temp) ; write out the original azimuth err = write_floats(outdir+'a_azm0',azm) ; write out field orientation, in degrees err = write_floats(outdir+'a_azm',azm) err = write_floats(outdir+'a_psi',psi) ;err = write_floats(outdir+'a__cct',cct) ; write the observed continuum instead of inversion value err = write_floats(outdir+'a__cct',occt) err = write_floats(outdir+'a_fld',fld) err = write_floats(outdir+'a_1azm',az0) err = write_floats(outdir+'a_2azm',az1) err = write_floats(outdir+'a_1incl',el0) err = write_floats(outdir+'a_2incl',el1) err = write_floats(outdir+'a__longitude',slong) err = write_floats(outdir+'a__latitude',slat) err = write_floats(outdir+'a__rgtasn',xx) err = write_floats(outdir+'a__dclntn',yy) err = write_floats(outdir+'a_dop',dop) err = write_floats(outdir+'a_dmp',dmp) err = write_floats(outdir+'a_eta0',eta0) err = write_floats(outdir+'a_b1mu',b1mu) err = write_floats(outdir+'a_bzero',bzero) err = write_floats(outdir+'a_alpha',alpha) ;err = write_floats(outdir+'a_delta',delta) ; write some auxiliary indices ;err = write_floats(outdir+'a__pvpv',signed_v) ;err = write_floats(outdir+'a__pip',netpol) err = write_floats(outdir+'a_chi',chisq) ;err = write_floats(outdir+'a_chi_ii',chisq_i) ;err = write_floats(outdir+'a_chi_qq',chisq_q) ;err = write_floats(outdir+'a_chi_uu',chisq_u) ;err = write_floats(outdir+'a_chi_vv',chisq_v) end