; this program computes the Doppler velocity and magnetogram ; using the same MDI-like algorithm as HMI ; from Stokes I simulated data from Kok-Leng Yeo ;------------------------------------------------------------------ PRO observables_KokLeng,iii MURAM=1 ; if MURAM=1 WE USE THE MURAM LINE PROFILE, OTHERWISE WE USE CO5BOLD sunspot=1 ; takes the Fe I line profile averaged over all sunspots' pixels ; READ THE HMI FILTER PROFILES AND Fe I LINE PROFILES ; MODIFY THIS PART ACCORDING TO THE ACTUAL NAMES YOU WILL USE RESTORE,'filter.sav' ; IDL save file that contains the 6 HMI filter profiles ; the code assumes that the filters ; are an array called filter with the ; following format: ; filter[wavelength,number] where ; number should be 6. The wavelengths ; are in Angstroms, relative to ; 6173.3433 A RESTORE,'average_profile.sav' ; IDL save file of the average Fe I line profile in a quiet region ; the code assumes that this file ; contains the line ; profile as an array called ; line_from_simulation, and the ; corresponding wavelengths as an ; array called ; wavelength_from_simulation. The ; wavelengths are in Angstroms, and 0 ; corresponds to the center of the Fe ; I line at rest (should be 6173.3433 A) ; PRODUCES THE HMI LOOK-UP TABLES N=6 ; 6 filters angle=DBLARR(N) angle[0]=+2.5d0 ; tuning positions in the order I0, I1, I2, I3, I4, and I5 (FROM RED WING TO BLUE WING) angle[1]=+1.5d0 ; angle[2]=+0.5d0 ; angle[3]=-0.5d0 ; angle[4]=-1.5d0 ; angle[5]=-2.5d0 ; FSR = 0.1689d0; FSR OF NB MICHELSON IN ANGSTROMS dtune = FSR/2.5d0 tune = angle*dtune angle = angle*2.d0*!dpi/DOUBLE(N) cosi = COS(angle) sini = SIN(angle) cos2i = COS(2.d0*angle) sin2i = SIN(2.d0*angle) dlamdv = 2.059205672212074294d-5 dvtune = dtune/dlamdv pv1 = dvtune*DOUBLE(N-1) pv2 = pv1/2.d0 magnetic = 1.d0/(2.d0*4.67d-5*0.000061733433d0*2.5d0*299792458.d0); //Lande factor=2.5 for Fe I line ntest = 821 ; number of test velocities for the look-up tables dvtest = 24.d0 ; sample rate for the test velocities in m/s vtest = (DINDGEN(ntest)-DOUBLE(ntest-1)/2.d0)*dvtest ; test velocities dlam = dvtest*dlamdv nlam = 32001 lam = (DINDGEN(nlam)-DOUBLE(nlam-1)/2.d0)*dlam lam0 = 6173.3433d0 ; YOU NEED TO INTERPOLATE THE FE I LINE PROFILE YOU HAVE ONTO THE lam ; GRID OF WAVELENGTHS. MODIFY THIS PART ACCORDINGLY. lineprofile = INTERPOL(line_from_simulation,wavelength_from_simulation,lam,/QUADRATIC) ; !!!!!! WARNING, I NARROWED THE LINE !!!!! a=WHERE(lam LT MIN(wavelength_from_simulation) OR lam GT MAX(wavelength_from_simulation)) lineprofile[a]=1.d0 a=WHERE(lineprofile eq MIN(lineprofile)) lam=lam-lam[a[0]] ; center the line so that 0 corresponds to the minimum of the line ; NOW YOU NEED TO INTERPOLATE THE FILTER PROFILES ONTO THE SAME lam GRID filters=DBLARR(nlam,N) FOR i=0,N-1 DO filters[*,i]=INTERPOL(REFORM(filter[*,i]),lamtemp,lam,/QUADRATIC) a=WHERE(lam LT MIN(lamtemp) OR lam GT MAX(lamtemp)) filters[a,*] = 0.d0 maxshiftlam=ROUND(vtest[ntest-1]*dlamdv/dlam) inten=DBLARR(N) vel1 =DBLARR(ntest) vel2 =vel1 ;GOTO,tables_computed ; if look-up tables have already been computed ;(UNCOMMENT THIS LINE, SO THAT YOU DON'T RECOMPUTE THE TABLES EVERY ;TIMEYOU RUN THE CODE) ; PLOT the filters and line profile for a quick look ; to make sure they both look fine plot,lam,lineprofile,xrange=[-0.3,0.3],xst=1,thick=2 oplot,lam,filters[*,0]/max(filters) oplot,lam,filters[*,1]/max(filters) oplot,lam,filters[*,2]/max(filters) oplot,lam,filters[*,3]/max(filters) oplot,lam,filters[*,4]/max(filters) oplot,lam,filters[*,5]/max(filters) ; loop over the test velocities FOR i=0,ntest-1 DO BEGIN shiftlam=ROUND(vtest[i]*dlamdv/dlam) ;SIGN CONVENTION: POSITIVE VELOCITIES CORRESPOND TO REDSHIFT PRINT,vtest[i] oplot,lam,shift(lineprofile,-shiftlam),col=180,thick=2 FOR j=0,N-1 DO BEGIN inten[j] = 0.d0 FOR k=maxshiftlam,nlam-maxshiftlam-1 DO inten[j] = inten[j]+filters[k,j]*lineprofile[k-shiftlam] ENDFOR f1c=0.d0 f1s=0.d0 f2c=0.d0 f2s=0.d0 FOR j=0,N-1 DO BEGIN f1c = f1c + cosi[j] *inten[j] f1s = f1s + sini[j] *inten[j] f2c = f2c + cos2i[j] *inten[j] f2s = f2s + sin2i[j] *inten[j] ENDFOR vel1[i] = atan(-f1s,-f1c)*pv1/2.d0/!dpi ; LOOK-UP TABLE FOR 1ST FOURIER VELOCITY vel2[i] = atan(-f2s,-f2c)*pv2/2.d0/!dpi ; LOOK-UP TABLE FOR SECOND FOURIER VELOCITY vel2[i] = ((vel2[i]-vel1[i]+10.5*pv2) MOD pv2)-pv2/2.0+vel1[i] ENDFOR SAVE,vtest,vel1,vel2,file="lookup.bin" tables_computed: RESTORE,"lookup.bin" ; READ THE SIMULATED FILTERGRAMS FROM KOK-LENG ; modify this section according to the name of your files, and their format HMI = READFITS(filename) ; Here I assume HMI is Stokes I ; if your simulation had a magnetic ; field, i.e. if I+V is different from ; I-V, then you need to read ; separately I+V and I-V, and then ; apply the MDI-like algorithm twice ; (once on I+V and once on I-V) and ; compute the average result for the ; Doppler velocity, linewidth, ; linedepth, continuum. The field ; strength is the difference in I+V ; velocity and I-V velocity multiplied ; by the variable called magnetic (I think...) savename = 'result.sav' ; dimensions of your 2D map (images). ; here i assume that your data are in a 2D array called HMI[x,y] ; you have to change the nx and ny values according to your maps nx = 600 ny = 600 Doppler = DBLARR(nx,ny) ; Doppler velocity Doppler2 = Doppler ; Doppler obtained from the 2nd Fourier coeffs. (not used in actual HMI data!) widthg = Doppler ; linewidth Idg = Doppler ; linedepth I0g = Doppler ; continuum intensity ; APPLY THE MDI-LIKE ALGORITHM N=double(N) period = (N-1.d0)*dtune FWHM=100.67102d0 ; reference FWHM of the Fe I line at disk center (might need to use another value if your reference line profile has a different FWHM FOR i=0,nx-1 DO BEGIN PRINT,i FOR j=0,ny-1 DO BEGIN f1c = TOTAL(cosi*REFORM(HMI[i,j,*])) ; 1st Fourier coefficient,cosine f1s = TOTAL(sini*REFORM(HMI[i,j,*])) ; 1st Fourier coefficient, sine f2c = TOTAL(cos2i*REFORM(HMI[i,j,*])); 2nd Fourier coefficient, cosine f2s = TOTAL(sin2i*REFORM(HMI[i,j,*])); 2nd Fourier coefficient, sine velocity = atan(-f1s,-f1c)*pv1/2.d0/!dpi velocity2= atan(-f2s,-f2c)*pv2/2.d0/!dpi velocity2= ((velocity2-velocity+10.5*pv2) MOD pv2)-pv2/2.0+velocity Doppler[i,j] = INTERPOL(vtest,vel1,velocity) ; WARNING: you might have to add a minus sign if your velocity convention is the opposite of HMI Doppler2[i,j]= INTERPOL(vtest,vel2,velocity2) temp = period/!dpi*sqrt(1.d0/6.d0*alog((f1c*f1c+f1s*f1s)/(f2c*f2c+f2s*f2s))) ; widthg[i,j] = (temp*2.d0*sqrt(alog(2.d0)))*1000.d0 ; we want the FWHM not the sigma of the Gaussian, in milliAngstoms Idg[i,j] = period/2.d0*sqrt(f1c*f1c+f1s*f1s)/sqrt(!dpi)/FWHM*exp(!dpi*!dpi*FWHM*FWHM/period/period) ; does not use the second Fourier coefficient temp3=Doppler[i,j]/dv I0g[i,j]=TOTAL(REFORM(HMI[i,j,*]))/double(N) for ii=0,N-1 DO I0g[i,j] += (Idg[i,j]/double(N)*exp(-(tune[ii]-temp3)^2.d0/FWHM/FWHM)) ENDFOR ENDFOR !p.multi=[0,2,2] TVIM,Doppler,/scale TVIM,field,/scale TVIM,Idg TVIM,I0g SAVE,doppler,doppler2,widthg,I0g,file=savename END