FUNCTION eis_density, fit1, fit2, ratio, line1=line1, line2=line2, $ isothermal=isothermal, no_yshift=no_yshift, quiet=quiet, $ yshift_wvl=yshift_wvl, $ ioneqfile=ioneqfile, abundfile=abundfile, $ calib=calib ;+ ; NAME ; ; EIS_DENSITY ; ; PROJECT ; ; Hinode/EIS ; ; EXPLANATION ; ; This routine takes as input two fit structures produced by ; eis_auto_fit and assumes that they correspond to the denominator ; and numerator lines of a density diagnostic. With the additional ; input of a theoretical line ratio curve, the routine computes ; the density at each pixel in the raster. ; ; A column depth is also calculated for each pixel using the ; emission measure method. The alternative isothermal method is ; also available as an option. ; ; For more details about densities, column depths and filling ; factors please consult EIS Software Note #15 which is available ; within Solarsoft at: ; ; $SSW/hinode/eis/doc/eis_notes/15_DENSITY/eis_swnote_15.pdf ; ; INPUTS ; ; FIT1 A structure in the format output by eis_auto_fit, ; representing the denominator line of the density ; sensitive line ratio. ; ; FIT2 A structure in the format output by eis_auto_fit, ; representing the numerator line of the density ; sensitive line ratio. ; ; RATIO A structure in the format output by ; eis_chianti_dens_ratio representing the theoretical ratio ; of the two lines as a function of density. ; ; OPTIONAL INPUTS ; ; LINE1 If the FIT1 structure contains fits for multiple lines, ; then LINE1 specifies the index of the line to be used. ; ; LINE2 If the FIT2 structure contains fits for multiple lines, ; then LINE2 specifies the index of the line to be used. ; ; NO_YSHIFT By default the routine automatically shifts the ; intensity array of the numerator line in the Y ; direction by the pixel amount given by ; EIS_CCD_OFFSET. Setting this keyword prevents the ; shift from being performed. ; ; YSHIFT_WVL This allows a specific wavelength to be specified ; through which the intensity arrays are aligned with ; each other. For example, if you were studying the Fe ; XIV 264/274 ratio, but wanted to align the output ; density map with Fe XII ; 195, then you would specify yshift_wvl=195.12. ; ; IONEQFILE The name of an ion balance file to use when computing ; the column depth. The default is to use the standard ; ion balance file in CHIANTI (!ioneq_file). ; ; ABUNDFILE The name of an element abundance file to use when ; computing the column depth. The default is to use the ; 'sun_coronal_ext.abund' file in CHIANTI. If ; this doesn't exist then the !abund_file is ; used (a photospheric abundance file). ; ; CALIB This allows the EIS calibration to be modified, changing ; the output intensity and error arrays. One of the ; following values should be set: ; 0 - no effect ; 1 - Del Zanna (2013) calibration is implemented. ; 2 - Del Zanna (2013) calibration, and the ; /correct_sensitivity correction is "undone". ; 3 - Warren et al. (2014) calibration is implemented. ; 4 - Warren et al. (2014) calibration, and the ; /correct_sensitivity correction is "undone". ; ; WARNING: before using this keyword, the user should know ; what calibration options were applied in the original ; call to eis_prep. See the EIS wiki for more details ; (http://solarb.mssl.ucl.ac.uk:8080/eiswiki/Wiki.jsp?page=EISCalibration) ; ; KEYWORDS ; ; ISOTHERMAL Two approximations are available when computing the ; column depth. The default is the emission measure ; approach. Setting /ISOTHERMAL will use the isothermal ; approach. A document is available on the EIS wiki ; describing the difference. ; ; QUIET If set, then information messaages will not be ; printed to the screen. ; ; OUTPUTS ; ; Returns a structure with the following tags: ; ; .ldens Logarithm (base 10) of the electron number density ; (cm^-3). ; .dens Electron number density (cm^-3). ; .lo The log of the lower limit to the density (cm^-3). ; .hi The log of the upper limit to the density (cm^-3). ; .coldepth The column depth of the plasma (cm). ; .coldepth_err Error on column depth (cm). ; .rat Value of the observed emission line ratio. ; .arcsec Scalar giving factor to convert cm to arcsec (for ; column depth). ; .abundfile Name of abundance file used in calculation. ; .ioneqfile Name of ion balance file used in calculation. ; .missing Value for missing data used in the arrays. ; .calib Stores the option used for the CALIB keyword. ; ; Each of the above quantities has the same dimensions as the ; intensity arrays of the input fit structures. ; ; Note that the .arcsec tag is used to convert the column depth ; values to arcsec units: ; ; coldepth_arcsec=dens.coldepth*dens.arcsec ; ; PROGRAMMING NOTES: ; The routine takes the intensity and intensity-error from the ; .int and .interr tags of the fitdata structures. It does not ; make use of the parameters in the .aa tag. So if you want to ; correct for a blend prior to calling eis_density, you can do it ; by modifying the .int and .interr tags only. ; ; CALLS ; ; ERRORS, EIS_GET_FITDATA, READ_IONEQ, READ_ABUND, CONVERTNAME, ; INTEGRAL_CALC ; ; HISTORY ; ; Beta 1, 18-Mar-2010, Peter Young ; Beta 2, 11-Aug-2010, Peter Young ; added check on muti-Gauss fits; now prints the theoretical ; ratio lines and observed lines so user can manually check that ; correct lines are being used. ; Ver.1, 5-May-2011, Peter Young ; routine now adjusts intensity array of numerator for CCD ; Y-offset; introduced /NO_YSHIFT in case user does not want to ; use this option. Added /QUIET keyword. Implemented missing ; data value for all arrays, and added 'missing' as a tag to the ; output structure. ; Ver.2, 6-May-2011, Peter Young ; Now calculates the column depth error and stores it in the ; new .coldepth_err tag. Also added .arcsec tag. ; Ver.3, 31-May-2011, Peter Young ; For emission measure method of calculating column depth, the ; density dependence of the integral is now included. ; Ver.4 25-Mar-2013, Peter Young ; The routine now makes use of the new 'quick' tag in the RATIO ; input in order to speed up the call to INTEGRAL_CALC. This ; makes a significant performance improvement for the coronal ; iron ions. ; Ver.5, 24-Jan-2014, Peter Young ; Added yshift_wvl= keyword. ; Ver.6, 16-Jul-2014, Peter Young ; Added abundfile= and ioneqfile= keywords; also added these as ; tags to output to structure; yshift_wvl has also been added to ; output. ; Ver.7, 30-Jun-2015, Peter Young ; Added CALIB= keyword for modifying the intensity calibration, ; and added calib tag to output. ; Ver.8, 24-Oct-2015, Peter Young ; Added information to header; no change to code. ; Ver.9, 5-Dec-2019, Peter Young ; Updated default coronal abundance file to latest CHIANTI ; file. ; Ver.10, 28-oct-2024, Peter Young ; Now handles the case where the measured ratio is outside of the ; density sensitivity range of the theoretical ratio. ;- IF n_params() LT 3 THEN BEGIN print,'Use: IDL> dens = eis_density( fit1, fit2, ratio [, line1=, line2= /isothermal,' print,' yshift_wvl=, /quiet, /no_yshift, abund_file=, ioneq_file=,' print,' calib= ])' print,' fit1 - auto_fit structure for denominator line' print,' fit2 - auto_fit structure for numerator line' print,' ratio - structure returned by eis_chianti_dens_ratio.pro' return,-1 ENDIF IF fit1.nx NE fit2.nx OR fit1.ny NE fit2.ny THEN BEGIN print,'%EIS_DENSITY: the array sizes of FIT1 and FIT2 are incompatible. Please check your inputs.' print,' Returning...' return,-1 ENDIF nx=fit1.nx ny=fit2.ny IF fit1.ngauss GT 1 THEN BEGIN IF n_elements(line1) EQ 0 THEN BEGIN print,'Fit structure 1 contains '+trim(fit1.ngauss)+' Gaussians. Please specify which Gaussian' print,'to use by specifying the LINE1 input.' print,'Returning...' return,-1 ENDIF ENDIF ELSE BEGIN line1=0 ENDELSE IF fit2.ngauss GT 1 THEN BEGIN IF n_elements(line2) EQ 0 THEN BEGIN print,'Fit structure 2 contains '+trim(fit2.ngauss)+' Gaussians. Please specify which Gaussian' print,'to use by specifying the LINE2 input.' print,'Returning...' return,-1 ENDIF ENDIF ELSE BEGIN line2=0 ENDELSE IF n_elements(calib) EQ 0 THEN calib=0 ; ; The YSHIFT_WVL keyword shifts INT2 in Y in order to ; spatially match INT1. Sub-pixel shifts are allowed for. ; IF n_elements(yshift_wvl) NE 0 THEN BEGIN int1=eis_get_fitdata(fit1,/int,err=int1_err,line=line1,yshift_wvl=yshift_wvl,calib=calib) int2=eis_get_fitdata(fit2,/int,err=int2_err,line=line2,yshift_wvl=yshift_wvl,calib=calib) ENDIF ELSE BEGIN IF keyword_set(no_yshift) THEN yshift_wvl=0. ELSE yshift_wvl=fit1.refwvl[line1] int1=eis_get_fitdata(fit1,/int,err=int1_err,line=line1,calib=calib) int2=eis_get_fitdata(fit2,/int,err=int2_err,line=line2,yshift_wvl=yshift_wvl,calib=calib) ENDELSE wvl1=fit1.refwvl[line1] wvl2=fit2.refwvl[line2] IF NOT keyword_set(quiet) THEN BEGIN print,'' print,'Theoretical ratio is: ',trim(ratio.num_wvl)+'/'+trim(ratio.den_wvl) print,'Observed ratio is: ',trim(string(format='(f12.3)',wvl2))+'/'+trim(string(format='(f12.3)',wvl1)) print,'' ENDIF IF n_elements(ioneqfile) NE 0 THEN BEGIN chck=file_search(ioneqfile) IF chck[0] EQ '' THEN BEGIN print,'% EIS_DENSITY: specified ion balance file not found. Please check your inputs.' print,' Returning...' return,-1 ENDIF ENDIF ELSE BEGIN ioneqfile=!ioneq_file ENDELSE ; ; Compute ratio and error ; errors, int2, int2_err, int1, int1_err, rat, err missing=fit1.missing i_miss=where(int2 EQ missing OR int1 EQ missing,n_miss) IF n_miss NE 0 THEN rat[i_miss]=missing rat_theory=ratio.ratio dens=ratio.dens em_sum=ratio.em_sum ; ; ratio_type=1 if the ratio increases with density, and 0 if it decreases. ; IF rat_theory[-1] GT rat_theory[0] THEN ratio_type=1b ELSE ratio_type=0b ; ; Interpolate onto 0.001 dex intervals ; nd=(max(dens)-min(dens))*1000.+1 densi=findgen(nd)/1000.+dens[0] y2=spl_init(dens,rat_theory) rati=spl_interp(dens,rat_theory,y2,densi) ; y2=spl_init(dens,em_sum) em_sumi=spl_interp(dens,em_sum,y2,densi) ; ; Re-define density values ; rat_theory=rati dens=densi em_sum=em_sumi outdens=make_array(nx,ny,/float,value=missing) hidens=outdens lodens=outdens ;hidens=fltarr(nx,ny) ;lodens=fltarr(nx,ny) emiss=fltarr(nx,ny) ; ; Go through each pixel in raster and compute density and the upper ; and lower limits. The emissivity sum at each location is also stored ; as this used to compute the isothermal column depth. ; ; 28-Oct-2024: I now check if the ratio is above or below the ; sensitivity range. If yes, then density is set to missing. If the ; upper or lower limit on the ratio lies within the sensitivity range ; then lodens or hidens will be set ; FOR i=0,nx-1 DO BEGIN FOR k=0,ny-1 DO BEGIN CASE 1 OF ; rat[i,k] GE max(rat_theory): BEGIN rat_err=rat[i,k]-err[i,k] IF rat_err LT max(rat_theory) THEN BEGIN getmin=min(abs(rat_theory-(rat_err)),j) IF ratio_type EQ 1 THEN BEGIN lodens[i,k]=dens[j] ENDIF ELSE BEGIN hidens[i,k]=dens[j] ENDELSE ENDIF END ; rat[i,k] LE min(rat_theory): BEGIN rat_err=rat[i,k]+err[i,k] IF rat_err GT min(rat_theory) THEN BEGIN getmin=min(abs(rat_theory-(rat_err)),j) IF ratio_type EQ 1 THEN BEGIN hidens[i,k]=dens[j] ENDIF ELSE BEGIN lodens[i,k]=dens[j] ENDELSE ENDIF END ; ELSE: BEGIN getmin=min(abs(rat_theory-rat[i,k]),j) outdens[i,k]=dens[j] ; ; Store emissivity sum for pixel ; emiss[i,k]=em_sum[j] ; getmin=min(abs(rat_theory-(rat[i,k]+err[i,k])),j) hidens[i,k]=dens[j] ; getmin=min(abs(rat_theory-(rat[i,k]-err[i,k])),j) lodens[i,k]=dens[j] END ENDCASE ENDFOR ENDFOR ;; ;; THE FOLLOWING SECTION OF CODE IS FOR COMPUTING THE COLUMN DEPTH ;; -------------------------------------------------------------- convertname,ratio.ion,iz,ion ; ; This indicates if the quick method for integral_calc should be used. ; IF tag_exist(ratio,'quick') EQ 1 THEN BEGIN quick=ratio.quick ENDIF ELSE BEGIN quick=0 ENDELSE ; ; The following loads up the element abundance. If abundfile= has been ; specified, then is used. Otherwise the software tries to use the ; coronal abundance file in CHIANTI. If this is not found for some ; reason, then the routine defaults to !abund_file which is a standard ; photospheric abundance file. ; IF n_elements(abundfile) NE 0 THEN BEGIN chck=file_search(abundfile) IF chck[0] EQ '' THEN BEGIN print,'% EIS_DENSITY: specified abundance file not found. Please check your inputs.' print,' Returning...' return,-1 ENDIF ENDIF ELSE BEGIN ab_dir=file_dirname(!abund_file) abundfile=concat_dir(ab_dir,'sun_coronal_1992_feldman_ext.abund') chck=file_search(abundfile) IF chck EQ '' THEN abundfile=!abund_file ENDELSE read_abund,abundfile,ab,ref ab=ab[iz-1] IF NOT keyword_set(quiet) THEN BEGIN print,'Column depth calculation' print,' Abundance file: '+trim(abundfile) print,' Ion balance file: '+trim(ioneqfile) print,'' ENDIF IF NOT keyword_set(isothermal) THEN BEGIN ; ; ---- Compute column depth using emission measure method ---- ; ; The density enters into the emission measure integral and so it is ; necessary to compute the integral for each pixel individually. This ; is very time consuming so the solution is to compute the integral ; for a range of densities between the MIN(DENS) and MAX(DENS), spaced ; at 0.5 dex intervals. The returned function will be a smoothly ; varying function of density and so I interpolate onto a very fine density ; scale. I then go through each pixel of the OUTDENS array and take ; the integral value that corresponds to the density in that pixel. ; This speeds up the process greatly while keeping the necessary ; accuracy. ; ; 25-Mar-2013: since CHIANTI 7.1, the call to integral_calc has become ; extremely slow due to the very large models for the coronal iron ; ions. To fix this I've introduced a quick method in ; integral_calc (see /quick keyword). ; ; Compute integral for sparse density intervals ; n_int=round((max(dens)-min(dens))/0.5)+1 dens_int=findgen(n_int)*0.5+min(dens) int_calc=dblarr(n_int) ; IF quick EQ 1 AND NOT keyword_set(quiet) THEN BEGIN print,'' print,' **Warning: this calculation may take several minutes.**' print,'' ENDIF ; FOR i=0,n_int-1 DO BEGIN integral_calc,iz,ion,dens=dens_int[i],index=[ratio.num_ind,ratio.den_ind], $ outstr=outstr, quick=quick, /quiet int_calc[i]=outstr.pidec ENDFOR emiss_out=0 ; ; Interpolate integral onto much finer density scale. ; nd=round((max(dens)-min(dens))*1000.+1) densi=findgen(nd)/1000.+dens[0] y2=spl_init(dens_int,int_calc) int_calc_i=spl_interp(dens_int,int_calc,y2,densi) ; ; Compute integral array pixel-by-pixel using nearest density ; ic_array=fltarr(nx,ny)+missing FOR i=0,nx-1 DO BEGIN FOR j=0,ny-1 DO BEGIN IF outdens[i,j] NE missing THEN BEGIN getmin=min(abs(outdens[i,j]-densi),k) ic_array[i,j]=int_calc_i[k] ENDIF ENDFOR ENDFOR ; ; Compute column depth ; coldepth=(int1+int2)*ic_array/ab/(10.^outdens)^2 ; ; Set missing data ; IF n_miss NE 0 THEN coldepth[i_miss]=missing ENDIF ELSE BEGIN ;------- Compute column depth using isothermal approximation ------ ; ; This approach assumes the plasma is isothermal at the T_max of the ; ion. ; ; ; Determine the ionization fraction (f_of_t) for the T_max of the ion ; read_ioneq,ioneqfile,tt,ii,ref ii=reform(ii[*,iz-1,ion-1]) f_OF_t=max(ii) ; alpha=4.0*!pi/0.83/ab/f_OF_t /emiss ; coldepth=alpha*(int1+int2)/10.^outdens ; IF n_miss NE 0 THEN coldepth[i_miss]=missing ENDELSE ; ; Compute column depth error (see EIS Software Note #15) ; sig_i2=int1_err^2+int2_err^2 a=-2.0*alog(10.) sig_n= ( ( hidens-outdens ) + ( outdens-lodens ) )/2. coldepth_err=sqrt(sig_i2*coldepth^2/(int1+int2)^2 + sig_n^2 * a^2 * coldepth^2) IF n_miss NE 0 THEN coldepth_err[i_miss]=missing dens=10.^outdens IF n_miss NE 0 THEN BEGIN dens[i_miss]=missing lodens[i_miss]=missing hidens[i_miss]=missing ENDIF dens={ldens: outdens, $ dens: dens, $ lo: lodens, $ hi: hidens, $ ratio: rat, $ coldepth: coldepth, $ coldepth_err: coldepth_err, $ arcsec: 1./7.2527e7, $ abundfile: abundfile, $ ioneqfile: ioneqfile, $ yshift_wvl: yshift_wvl, $ calib: calib, $ missing: missing} return,dens END