;+ ; Project : SOHO - CDS ; ; Name : DMM_CALC_SPEC ; ; Purpose : Heart of Dere et al (CHIANTI) synthetic spectrum calculation. ; ; Explanation : ; ; Use : Mainly from within procedure CHIANTI_SS, but could stand alone ; as ; IDL> dmm_calc_spec, input_str, output_str ; ; Inputs : input_str - structure with details of the setup. Has tags: ; isg - the spectrograph 0-GIS, 1-NIS) ; ichannel - which spectral band (0-3 or 0-1) ; slith - slit height in arcseconds ; slitw - slit width in arcseconds ; slit_id - character string slit description ; abund_name - full path and name of abundance file ; dem_name - full path and name of dem file ; ioneq_name - full path and name of ioneq file ; pressure - pressure to use in the calculation ; wvlmin - minimum wavelength ; wvlmax - maximum wavelength ; ; Opt. Inputs : None ; ; Outputs : None ; ; Opt. Outputs: output_str - structure giving calculated spectrum etc. Tags: ; ; lambda - array of wavelengths of spectrum ; spectrum - array of intensities ; list_wvl - wavelengths of lines used ; list_ident - identifications of lines used ; nlines - number of spectral lines used. ; ; Keywords : None ; ; Calls : setup_ion ; ; Common : None ; ; Restrictions: None ; ; Side effects: None ; ; Category : Spectrum ; ; Prev. Hist. : Based on cds_spectrum by K. Dere (NRL) ; ; Written : CDS version, linked to dmm_ss by C D Pike, RAL, 14-Jul-95 ; ; Modified : Pick up wavelength range from input structure. CDP, 19-Jul-95 ; Changes to calculation. KD, 14-Feb-96 ; Data file cds_nis_wvl.dat changed to reflect actual wavelength ; coverage. CDP, 11-Jun-96 ; Move auxiliary .dat files to !PATH. CDP, 12-Sep-96 ; Rework to cut out all instrumental effects - apart ; from broadening. CDP, 27-Sep-96 ; ; Modified : Version 8, 23-Sep-97 ; Changed to take into account the new DEM file format ; (log(DEM) values instead of log(T*DEM) values ) ; Giulio Del Zanna, UCLAN, 27-Apr-97 ; ; Modified : Version 9, 5-Oct-97 ; Changed few minor things, mainly how the density is ; calculated. To uniform with Synthetic. ; Differences of about 5% ,compared with ; the previous version, are to be found. ; Also changed the UNITS of the output. Now in ergs/cm2/s/str. ; GDZ, UCLAN, 5-Oct-97 ; V.10, Update element list, CDP, 18-Jun-99 ; V. 11, Giulio Del Zanna (DAMTP), 30-Oct-2000. ; Rewritten completely the routine, to make it compatible with ; CHIANTI v.3 ; V. 12, reduced the array size to 10,000 lines max. ; GDZ 20-Nov-2000 ; ; Version : Version 12, 20-Nov-2000 ; ;- pro dmm_calc_spec, setup_struct, return_struct ; ; setup comes in as a structure for convenience so convert to individual ; variables to match original program ; isg = setup_struct.isg ichannel = setup_struct.ichannel slith = setup_struct.slith slitw = setup_struct.slitw slit_id = setup_struct.slit_id abund_name = setup_struct.abund_name dem_name = setup_struct.dem_name ioneq_name = setup_struct.ioneq_name pressure = setup_struct.pressure ; ; override standard settings with user choice ; wmin=setup_struct.wvlmin wmax=setup_struct.wvlmax ; ;derived from synthetic.pro ; ; wvlmin=short wavelength limit in Angstroms ; wvlmax=long wavelength limit in Angstroms ; pressure= electron pressure (Ne x T) cm-3 K ; ; ; ; commons with setup_ion: common elvlc,l1a,term,conf,ss,ll,jj,ecm,eryd,ecmth,erydth,eref common wgfa, wvl,gf,a_value common upsilon,t_type,deu,c_ups,splups ;common elements,abund,abund_ref,ioneq,ioneq_t,ioneq_ref spectrum = fltarr(long(wmax-wmin)*100.) lambda = findgen(long(wmax-wmin)*100.)*0.01 + wmin ; ; ; ; will be doing a convolution of the spectrum and then chopping off the ends ; read_abund,abund_name,abund,abund_ref ; read_ioneq,ioneq_name,ioneq_t,ioneq,ioneq_ref n_ioneq_t = n_elements(ioneq_t) dlnt=alog(10.^(ioneq_t(1)-ioneq_t(0))) ; read_dem,dem_name,dem_t,dem,dem_ref ; ; gdt=where((ioneq_t ge min(dem_t)) and (ioneq_t le max(dem_t))) dem_int1=10.^(spline(dem_t,dem,ioneq_t(gdt))) dem_int=fltarr(n_ioneq_t) ngt=n_elements(gdt) for igt=0,ngt-1 do begin dem_int(gdt(igt))=dem_int1(igt) endfor ; ; element=['H','He','Li','Be','B ','C ','N ','O ','F ','Ne','Na',$ 'Mg','Al','Si','P ','S ','Cl','Ar','K ','Ca','Sc','Ti',$ 'V ','Cr','Mn','Fe','Co','Ni','Cu','Zn'] ; ionstage=['I','II','III','IV','V','VI','VII','VIII','IX','X','XI','XII',$ 'XIII','XIV','XV','XVI','XVII','XVIII','XIX','XX','XXI',' XXII',$ 'XXIII','XXIV','XXV','XXVI','XXVII','XXVIII'] ; spd=['S','P','D','F','G','H','I','K'] jvalue=['0','1/2','1','3/2','2','5/2','3','7/2','4','9/2','5','11/2',$ '13/2','15/2'] ; maxwvl=10000 ; mast_z=intarr(maxwvl) mast_ion=intarr(maxwvl) mast_snote=strarr(maxwvl) mast_lvl=intarr(maxwvl) mast_wvl=fltarr(maxwvl) mast_A=fltarr(maxwvl) ; mast_term1=strarr(80,maxwvl) mast_desig1=strarr(80,maxwvl) mast_term2=strarr(80,maxwvl) mast_desig2=strarr(80,maxwvl) ; mast_abund=fltarr(maxwvl) mast_ioneq=fltarr(maxwvl) mast_tmax=fltarr(maxwvl) mast_pop=fltarr(maxwvl) mast_popazidem=fltarr(maxwvl) mast_spectrum=lonarr(maxwvl) ; ; mast_index=0 ; ; ; ; open the file that has the names of the ions ; defsysv,'!xuvtop', EXISTS = EXISTS IF NOT EXISTS THEN use_chianti defsysv,'!xuvtop', EXISTS = EXISTS IF NOT EXISTS THEN $ message, 'system variable !xuvtop must be set using e.g. use_chianti,directory ' ;xuvtop = getenv('CDS_SS_DERE') mname=!xuvtop+'/masterlist/masterlist.ions' read_masterlist,mname,list nlist=n_elements(list) ; ; main input and calculation loop ************** ; for ilist=0,nlist-1 do begin gname=list(ilist) convertname,gname,iz,ion ion2spectroscopic,gname,snote ; read in level information, wavelengths, gf and A values from .wgfa files ; setup_ion,gname,wmin,wmax,wvltst,lvl1,lvl2,wvl1,gf1,a_value1 ; ; does this ion have lines in the specified wavelength range? ; if (wvltst gt 0) then begin ;*********************************************** ; ; locname=strlowcase(gname) pos=strpos(locname,'_') l=strlen(pos) first=strmid(locname,0,pos) last=strmid(locname,pos+1,l-pos-1) ; if strpos(last,'d') ge 0 then dielectronic=1 else dielectronic=0 ; ntrans=n_elements(lvl1) nlvls=max([lvl1,lvl2]) ; ; calculate level populations ; this_ioneq=ioneq(*,iz-1,ion-1+dielectronic) itmax=where(this_ioneq eq max(this_ioneq)) ltmax=ioneq_t(itmax) tmax=10.^(ioneq_t(itmax)) tmax=tmax(0) ; necessary to make tmax a scalar ; edensity=pressure/tmax ; ; limit range of temperatures to those where dem is specified ; gioneq_t=where((this_ioneq gt 0.) and (dem_int gt 0.),ngt) ; popazidem=dblarr(nlvls) ; ; if ngt gt 0 then begin for ig=0,ngt-1 do begin it=gioneq_t(ig) tt=10.^ioneq_t(it) dt=tt*dlnt ; ; if keyword_set(density) then begin ; edensity=density ; endif else edensity=pressure/tt edensity=pressure/tt pop_solver,tt,edensity,pop pop = reform(pop) popazidem(0)=popazidem+pop(*)*this_ioneq(it)*dem_int(it)*dt/edensity endfor for itrans=0,ntrans-1 do begin l1=lvl1(itrans)-1 l2=lvl2(itrans)-1 ww=wvl1(itrans) ; ; can the level be excited? upssiz=size(splups) upssiz2=upssiz(2)-1 if(l2 le upssiz2) then begin maxups=max(splups(*,l2,0)) if(maxups gt 0.) then upstst=1 else upstst=0 endif else upstst=0 ; ; select only observed wvl's test1=(ww gt wmin) and (ww le wmax) test2=upstst gt 0 if test1 and test2 then begin mast_z(mast_index)=iz mast_ion(mast_index)=ion mast_snote(mast_index)=snote mast_lvl(mast_index)=lvl2(itrans) mast_wvl(mast_index)=wvl1(itrans) mast_A(mast_index)=a_value1(itrans) mast_abund(mast_index)=abund(iz-1) ; mast_tmax(mast_index)=ltmax mast_ioneq(mast_index)=this_ioneq(itmax) mast_popazidem(mast_index)=popazidem(l2) mast_spectrum(mast_index)=wvl1(itrans) ;(wvl1(itrans)+dlambda/2.-wmin)/dlambda ; ; get lower level designation term0=strtrim(term(l1),2) ; term2='' blank=strpos(term0,' ') while blank gt 0 do begin term2=term2+' '+strmid(term0,0,blank) term0=strmid(term0,blank,100) term0=strtrim(term0,2) blank=strpos(term0,' ') endwhile mast_term1(mast_index)=strlowcase(term2) ; jinteger=fix(2.*jj(l1)) jstring=jvalue(jinteger) ; spins=strtrim(string(ss(l1),'(i2)'),2) mast_desig1(mast_index)=spins+spd(ll(l1))+jstring ; ; ; get upper level designation term0=strtrim(term(l2),2) ; term2='' blank=strpos(term0,' ') ; print,term0 while blank gt 0 do begin term2=term2+' '+strmid(term0,0,blank) term0=strmid(term0,blank,100) term0=strtrim(term0,2) blank=strpos(term0,' ') endwhile mast_term2(mast_index)=strlowcase(term2) ; jinteger=fix(2.*jj(l2)) jstring=jvalue(jinteger) ; spins=strtrim(string(ss(l2),'(i2)'),2) mast_desig2(mast_index)=spins+spd(ll(l2))+jstring ; mast_index=mast_index+1 ENDIF ENDFOR ; over itrans (spectral lines) ; ENDIF ELSE BEGIN ; print,' ionization equilibrium eq 0 in DEM temperature interval' ; ngt gt 0 ENDELSE ENDIF ELSE BEGIN ; print,' no lines in specified wavelength interval' ; wvltest ENDELSE ; ENDFOR ; ; end of main input and calculation loop *************** ; ;print,' number of lines=',mast_index ; ; if mast_index ge 1 then begin mast_z=mast_z(0:mast_index-1) mast_ion=mast_ion(0:mast_index-1) mast_snote=mast_snote(0:mast_index-1) mast_lvl=mast_lvl(0:mast_index-1) mast_wvl=mast_wvl(0:mast_index-1) mast_A=mast_A(0:mast_index-1) mast_term1=mast_term1(0:mast_index-1) mast_desig1=mast_desig1(0:mast_index-1) mast_term2=mast_term2(0:mast_index-1) mast_desig2=mast_desig2(0:mast_index-1) ; mast_abund=mast_abund(0:mast_index-1) mast_tmax=mast_tmax(0:mast_index-1) mast_ioneq=mast_ioneq(0:mast_index-1) mast_popazidem=mast_popazidem(0:mast_index-1) mast_spectrum=mast_spectrum(0:mast_index-1) ; ; ; srt_index=sort(mast_wvl) mast_z(0)=mast_z(srt_index) mast_ion(0)=mast_ion(srt_index) mast_snote(0)=mast_snote(srt_index) mast_lvl(0)=mast_lvl(srt_index) mast_wvl(0)=mast_wvl(srt_index) mast_A=mast_A(srt_index) mast_term1=mast_term1(srt_index) mast_desig1=mast_desig1(srt_index) mast_term2=mast_term2(srt_index) mast_desig2=mast_desig2(srt_index) ; mast_abund=mast_abund(srt_index) mast_tmax=mast_tmax(srt_index) mast_ioneq=mast_ioneq(srt_index) mast_popazidem=mast_popazidem(srt_index) mast_spectrum=mast_spectrum(srt_index) ; list_ident=strarr(mast_index) list_wvl=fltarr(mast_index) ; nlines=mast_index ;print,' number of lines =',nlines ; ;for i=0,nlines-1 do begin ;print,mast_z(i),mast_ion(i),mast_lvl(i),mast_wvl(i),mast_A(i),mast_tmax(i), $ ; format='$(3i3,f12.3,2e10.3)' ;endfor ; lines=0 ; ; cspectrum = spectrum*0.0 for i=0,nlines-1 do begin ; ; elstage=element(mast_z(i)-1)+' '+ionstage(mast_ion(i)-1) wvls=strtrim(string(mast_wvl(i),'(f12.3)'),2) term1=strtrim(mast_term1(i),2) desig1=strtrim(mast_desig1(i),2) term2=strtrim(mast_term2(i),2) desig2=strtrim(mast_desig2(i),2) ; ; az=string(mast_abund(i),'(e9.2)') item1=string(mast_tmax(i),'(f4.1)') ; item2=string(mast_ioneq(i),'(e9.2)') ; item2=string(mast_pop(i),'(e9.2)') ; do intensity in photons/cm2/s/str ; intensity=mast_abund(i)*mast_A(i)*mast_popazidem(i)/(4.d*3.14159d) ; do intensity in ergs/cm2/s/str hc=6.626d-27*2.998d+10*1.d+8/(4.d*3.14159d*mast_wvl(i)) intensity=hc*mast_abund(i)*mast_A(i)*mast_popazidem(i) ; ; item2=string(intensity,'(e10.2)') ; list_wvl(i)=mast_wvl(i) list_ident(i)=mast_snote(i)+' '+term1+' '+desig1+' - '+term2+' '+desig2 list_ident(i)=strpad(list_ident(i),50,/after)+' Int='+item2+' Tmax='+item1 ; nn = round( (list_wvl(i)-wmin)*100) if nn lt n_elements(cspectrum) then cspectrum(nn) = cspectrum(nn)+intensity ; ; endfor ; ; tidy list_ident ; for i=0,n_elements(list_ident)-1 do begin nid1 =strpos(list_ident(i),'Int') text = strcompress(strmid(list_ident(i),0,45)) list_ident(i) = text + strmid(list_ident(i),45,150) n1 = strpos(list_ident(i),' ') n2 = strpos(list_ident(i),' ',n1+1) text = strpad(strmid(list_ident(i),0,n2),10,/after) list_ident(i) = text + strmid(list_ident(i),n2,150) nid2 =strpos(list_ident(i),'Int') text = strpad(strmid(list_ident(i),0,nid2),nid1+5,/after) list_ident(i) = text+strmid(list_ident(i),nid2,150) endfor ; ; fill out return structure ; return_struct = {lambda:lambda, spectrum:cspectrum, list_wvl:list_wvl,$ list_ident:list_ident, nlines:mast_index} ENDIF ELSE BEGIN return_struct = 0 endelse end