Документ взят из кэша поисковой машины. Адрес оригинального документа : http://star.arm.ac.uk/~jgd/outgoing/TEMP/gft_c_si_o_gd.pro Дата изменения: Tue Sep 8 11:16:14 2015 Дата индексирования: Sun Apr 10 06:27:22 2016 Кодировка: Поисковые слова: vallis

;----------------------------------------------------------------------
;
; PROJECT : ADAS - transient contribution functions studies for
; C IV, Si IV, O IV and O III
;
; NAME : gft_c_si_o
;
; PURPOSE : to calculate contribution functions as a function of
; electron temperature and electron density at equilibrium
; and as a function of electron temperature, electron
; density and integration time in the transient case
; for the following lines:
; C IV 1548.189 2s 2S1/2-2p 2P3/2
;
; Si IV 1393.757 3s 2S1/2-3p 2P3/2
;
; C III 977.02 2s2 1S0-2s 2p 1P1
; C III 1908.734 2s2 1S0-2s 2p 3P4
;
; O IV 1399.766 2s2 2p 2P1/2-2s 2p2 4P1/2
; O IV 1401.158 2s2 2p 2P3/2-2s 2p2 4P5/2
; O IV 1404.779 2s2 2p 2P3/2-2s 2p2 4P3/2
; O IV 1407.375 2s2 2p 2P3/2-2s 2p2 4P1/2
; O IV 787.71 2s2 2p 2P1/2-2s 2p2 2D3/2
; O IV 790.114 2s2 2p 2P3/2-2s 2p2 2D3/2
; O IV 790.201 2s2 2p 2P3/2-2s 2p2 2D5/2
;
; O III 832.929 2s2 2p2 3P0-2s 2p3 3D1
; O III 833.715 2s2 2p2 3P1-2s 2p3 3D1
; O III 833.749 2s2 2p2 3P1-2s 2p3 3D2
; O III 835.092 2s2 2p2 3P2-2s 2p3 3D2
; O III 835.289 2s2 2p2 3P2-2s 2p3 3D3
; O III 1666.150 2s2 2p2 3P4-2s 2p3 5S2
;
;
; INPUT : /c3 or keywords to choose the set of lines
; required /si3 or to be analysed. Once the ion is
; /o3 or selected, the routine will list
; /c2 or the corresponding lines available
; /o2 and ask for the respective index
;
;
; INPUT : density_index integer index correponding to
; optional one value in the density array
; if density_index is not provides, the routine
; will list the density values
; and ask to enter an index interactively
;
; /norec keyword if /norec is selected, the contribution
; due to recombination is ignored
;
; /exc keyword if /exc is selected, ic pecs for O3+ and
; Si3+ are used, inlcuding escitation
; contribution only,no projection
;
; tekin float initial temperature condition
; value for the transient (Kelvin)
; if tekin is not provides, the routine
; uses a default value= 30000. K
;
; densin float initial density condition
; value for the transient (cm-3)
; if densin is not provides, the routine
; uses a default value= 5.e10 cm-3
;
; OUTPUT : 1. plot output it is a .ps file which contains:
; - surface of equilibrium
; contribution functions and
; transient contribution
; functions vs electron temperature
; and integration time at a fixed
; value of electron density
; - plot of equilibrium (thick
; line) and transient contribution functions vs
; electron temperature at different times
; 2. save output it is a .save file which contains all
; variables: tekin float initial temperature
; densin float initial density
; tek array() final temperature
; 1st dim: num of temp.
; dens array() final density
; 1st dim: num of dens.
; time array() integration time
; 1st dim: num of time
; gnteq array(,) equilibrium contribution
; functions
; 1st dim: num of temp.
; 2nd dim: num of dens.
; gnt array(,,) transient contribution
; functions
; 1st dim: num of temp.
; 1st dim: num of dens.
; 1st dim: num of time
; 3. text output it is a .txt file which contains tekin,densin,tek,time,
; gnteq, gnt for a fixed value of final density
;
;
; CALLS : split_multiplet.pro, read_adf15.pro, run_adas405.pro, run_adas406.pro
;
; EXAMPLE : Examples of how to run:
;IDL> gft_c_si_o,/c3
;
;On the terminal will appear the following:
;
;% Compiled module: ADAS_VECTOR.
; *************************************
; INITIAL PLASMA CONDITIONS FOR TRANSIENT
; Electron temperature (K)= 30000.0
; Electron density (cm-3)= 1.00000e+10
; ************************************* ;
;
; Calculating contribution function for O IV 1401 Ang. or O IV 787 Ang.
;
; O IV 1401 2s2 2p 2P-2s 2p2 4P multiplet
; Index Line:
; 0 1399 2P1/2-4P1/2
; 1 1401 2P3/2-4P5/2
; 2 1404 2P3/2-4P3/2
; 3 1407 2P3/2-4P1/2
; O IV 787 2s2 2p 2P-2s 2p2 2D multiplet
; Index Line:
; 4 787 2P1/2-2D3/2
; 5 790 2P3/2-2D3/2
; 6 791 2P3/2-2D5/2
;Insert index for the line (e.g. 3): 4 <--this number is inserted by
; hand when requested
;% Compiled module: SPLIT_MULTIPLET.
;% Compiled module: READ_ADF15.
;% Compiled module: FILE_ACC.
;% Compiled module: ADAS_STRING_PAD.
;% Compiled module: NUMLINES.
;% Compiled module: RUN_ADAS405.
;% Compiled module: I4EIZ0.
;% Compiled module: XXDATE.
;% Compiled module: RUN_ADAS406.
; Index Electron density values (cm-3): <--this number is inserted by
; hand when requested
; 0 1.0000000e+10
; 1 1.5848932e+10
; 2 2.5118865e+10
; 3 3.9810719e+10
; 4 6.3095736e+10
; 5 1.0000000e+11
; 6 1.5848934e+11
; 7 2.5118863e+11
; 8 3.9810719e+11
; 9 6.3095727e+11
; 10 1.0000000e+12
;Insert index for density (e.g. 3): 5 <--this number is inserted by
; hand when requested
;
; PLOT OUTOOUT: o3_789.40_5_plot.ps
;
; SAVE OUTPUT: o3_789.40.save
;
; TEXT OUTPUT: o3_789.40_5_text.txt
;IDL>
;
;
;----------------------------------
; AUTHOR : Alessandra Giunta
;
; DATE : 14-04-2013
;
;
; MODIFIED : Alessandra Giunta
; 1.1
; - First version.
; 1.2 Added /exc keyword to allow the use of full ic pec
; but no full gcr for o3 and si3
;
;
; VERSION:
; 1.1 14-04-2013
; 1.2 09-12-2013
;
;
;----------------------------------------------------------------------
pro gft_c_si_o_gd,c3=c3,si3=si3,o3=o3,c2=c2,o2=o2,$
density_index=density_index,norec=norec,$
tekin=tekin,densin=densin,exc=exc


;define electron temperature, electron density and time arrays
tek=adas_vector(low=1.e4,high=1.e7,num=48)
te=tek/11605.
dens=adas_vector(low=1.e9,high=1.e13,num=11)
time=adas_vector(low=10.0,high=20.0,num=2,/linear)
;time=adas_vector(low=5.,high=100.,num=10,/linear)

nte=n_elements(te)
ndens=n_elements(dens)
ntime=n_elements(time)


tion=fltarr(nte)
trec=fltarr(nte)
tdion=fltarr(nte,ndens)
tdrec=fltarr(nte,ndens)
ion=fltarr(nte,ndens,ntime)
rec=fltarr(nte,ndens,ntime)

tioneq=fltarr(nte)
treceq=fltarr(nte)
ioneq=fltarr(nte,ndens)
receq=fltarr(nte,ndens)


gnt=dblarr(nte,ndens,ntime)
gnteq=dblarr(nte,ndens)

uid = 'adas'
if keyword_set(c3) or keyword_set(c2) then elem = 'c'
if keyword_set(si3) then elem = 'si'
if keyword_set(o3) or keyword_set(o2) then elem = 'o'
year = 96



;INITIAL CONDITIONS FOR TRANSIENT
;electron temperature
if not keyword_set(tekin) then tekin=3.e4
tein=tekin/11605.0
;electron density
if not keyword_set(densin) then densin=1.e8


print,' ************************************* '
print,' INITIAL PLASMA CONDITIONS FOR TRANSIENT '
print,' Electron temperature (K)= ',tekin
print,' Electron density (cm-3)= ',densin
print,' ************************************* '


;C IV 1548.1890 1s2 2s 2S1/2-1s2 2p 2P3/2
if keyword_set(c3) then begin
print,''
print,' Calculating contribution function for C IV 1548.1890 Ang.'
print,''
split_c3=split_multiplet(0.5, 0, 0.5, 1, /norm, j_low=j1, j_up=j2)

tt=total(split_c3)
t2=split_c3[1]/tt
t1=split_c3[0]/tt

;Use read_adf15.pro to read the PEC file for C+3
file='/home/adas/adas/adf15/pec96#c/pec96#c_pju#c3.dat'

block1=14
block2=55

;ionisation stages
is=3
ir=4
endif

;Si IV 1393.757 3s 2S1/2-3p 2P3/2
if keyword_set(si3) then begin
print,''
print,' Calculating contribution function for Si IV 1393.757 Ang.'
print,''

;Use read_adf15.pro to read the PEC file for Si+3
;---------------------------------------
if keyword_set(exc) then begin
file='pec_si3_lgy09.dat'
endif else begin
file='pec96#si_pju#si3.dat'
endelse
;---------------------------------------
if keyword_set(exc) then begin
t2=1.
block1=1
endif else begin
;---------------------------------------
split_si3=split_multiplet(0.5, 0, 0.5, 1, /norm, j_low=j1, j_up=j2)

tt=total(split_si3)
t2=split_si3[1]/tt
t1=split_si3[0]/tt

block1=12
block2=33
;---------------------------------------
endelse
;---------------------------------------

;ionisation stages
is=3
ir=4

endif

;O IV 1401 2s2 2p 2P-2s 2p2 4P and O IV 787 2s2 2p 2P-2s 2p2 2D multiplets
if keyword_set(o3) then begin
print,''
print,' Calculating contribution function for O IV 1401 Ang. or O IV 787 Ang.'
print,''


;PEC file for O+3

;---------NEW-09-12-13---------------
if keyword_set(exc) then begin
file='pec_o3_lgy12.dat'
endif else begin
file='peco3_780_1402.dat'
endelse
;------------------------------------

str=['1399 2P1/2-4P1/2','1401 2P3/2-4P5/2','1404 2P3/2-4P3/2','1407 2P3/2-4P1/2',$
'787 2P1/2-2D3/2','790 2P3/2-2D3/2','791 2P3/2-2D5/2']
nstr=n_elements(str)
print,' O IV 1401 2s2 2p 2P-2s 2p2 4P multiplet'
print,' Index Line:'
for l=0,3 do print,l,' ',str[l]
print,' O IV 787 2s2 2p 2P-2s 2p2 2D multiplet'
print,' Index Line:'
for l=4,nstr-1 do print,l,' ',str[l]
read,ll,prompt='Insert index for the line (e.g. 3): '
ll=fix(ll)
;---------------------------------------------
if keyword_set(exc) then begin
t2=1.
if ll eq 0 then begin
block1=6
endif
if ll eq 1 then begin
block1=4
endif
if ll eq 2 then begin
block1=5
endif
if ll eq 3 then begin
block1=7
endif
if ll eq 4 then begin
block1=1
endif
if ll eq 5 then begin
block1=2
endif
if ll eq 6 then begin
block1=3
endif

endif else begin
;---------------------------------------------
if ll ge 0 and ll le 3 then begin

block1=2
block2=52

split_o3=split_multiplet(0.5, 1, 1.5, 1, /norm, j_low=j1, j_up=j2)
tt=total(split_o3)
if ll eq 0 then begin
sp=where(j1 eq 0.5 and j2 eq 0.5)
sp=fix(sp[0])
t2=split_o3[sp]/tt
endif
if ll eq 1 then begin
sp=where(j1 eq 1.5 and j2 eq 2.5)
sp=fix(sp[0])
t2=split_o3[sp]/tt
endif
if ll eq 2 then begin
sp=where(j1 eq 1.5 and j2 eq 1.5)
sp=fix(sp[0])
t2=split_o3[sp]/tt
endif
if ll eq 3 then begin
sp=where(j1 eq 1.5 and j2 eq 0.5)
sp=fix(sp[0])
t2=split_o3[sp]/tt
endif
endif else begin

block1=1
block2=51

split_o3=split_multiplet(0.5, 1, 0.5, 2, /norm, j_low=j1, j_up=j2)
tt=total(split_o3)
if ll eq 4 then begin
sp=where(j1 eq 0.5 and j2 eq 1.5)
sp=fix(sp[0])
t2=split_o3[sp]/tt
endif
if ll eq 5 then begin
sp=where(j1 eq 1.5 and j2 eq 1.5)
sp=fix(sp[0])
t2=split_o3[sp]/tt
endif
if ll eq 6 then begin
sp=where(j1 eq 1.5 and j2 eq 2.5)
sp=fix(sp[0])
t2=split_o3[sp]/tt
endif
endelse

;---------------------------------------
endelse
;---------------------------------------


;ionisation stages
is=3
ir=4

endif

;C III 977 2s2 1S0-2s 2p 1P1 and C III 1908 2s2 1S0-2s 2p 3P4 multiplets
if keyword_set(c2) then begin
print,''
print,' Calculating contribution function for C III 977 Ang. or C III 1908 Ang.'
print,''


;PEC file for C+2
file='pecc2_970_1912.dat'

str=['977 1S0-1P1','1908 1S0-3P4']
nstr=n_elements(str)
print,' C III 977 2s2 1S0-2s 2p 1P1'
print,' Index Line:'
print,0,' ',str[0]
print,' C III 1908 2s2 1S0-2s 2p 3P4'
print,' Index Line:'
print,1,' ',str[1]
read,ll,prompt='Insert index for the line (e.g. 1): '
ll=fix(ll)


if ll eq 0 then begin

block1=1
block2=51

t2=1.

endif else begin

block1=3
block2=53

t2=1.

endelse

;ionisation stages
is=2
ir=3

endif


;O III 834 2s2 2p2 3P-2s 2p3 3D and O III 1666 2s2 2p2 3P-2s 2p3 5S multiplets
if keyword_set(o2) then begin
print,''
print,' Calculating contribution function for O III 834 Ang. or O III 1666 Ang.'
print,''


;PEC file for O+2
file='peco2_833_1700.dat'

str=['832.93 3P0-3D1','833.71 3P1-3D1','833.75 3P1-3D2',$
'835.09 3P2-3D2','835.29 3P2-3D3','1666.15 3P4-5S2']
nstr=n_elements(str)
print,' O III 834 2s2 2p2 3P-2s 2p3 3D multiplet'
print,' Index Line:'
for l=0,4 do print,l,' ',str[l]
print,' O III 1666 2s2 2p2 3P4-2s 2p3 5S multiplet'
print,' Index Line:'
for l=5,nstr-1 do print,l,' ',str[l]
read,ll,prompt='Insert index for the line (e.g. 3): '
ll=fix(ll)

if ll ge 0 and ll le 4 then begin

block1=1
block2=50

split_o2=split_multiplet(1., 1., 1., 2., /norm, j_low=j1, j_up=j2)
tt=total(split_o2)
if ll eq 0 then begin
sp=where(j1 eq 0. and j2 eq 1.)
sp=fix(sp[0])
t2=split_o2[sp]/tt
endif
if ll eq 1 then begin
sp=where(j1 eq 1. and j2 eq 1.)
sp=fix(sp[0])
t2=split_o2[sp]/tt
endif
if ll eq 2 then begin
sp=where(j1 eq 1. and j2 eq 2.)
sp=fix(sp[0])
t2=split_o2[sp]/tt
endif
if ll eq 3 then begin
sp=where(j1 eq 2. and j2 eq 2.)
sp=fix(sp[0])
t2=split_o2[sp]/tt
endif
if ll eq 4 then begin
sp=where(j1 eq 2. and j2 eq 3.)
sp=fix(sp[0])
t2=split_o2[sp]/tt
endif

endif else begin

block1=2
block2=51

t2=1.
endelse

;ionisation stages
is=2
ir=3

endif



read_adf15,file=file,block=block1,te=te,dens=dens, $
data=dats,wlngth=wlngth,/all

;pec excitation
datas=dats*t2

;------------------------------------
if not keyword_set(exc) then begin
;------------------------------------
read_adf15,file=file,block=block2,te=te,dens=dens, $
data=datr,wlngth=wlngth,/all

;pec recombination
datar=datr*t2
;------------------------------------
endif
;------------------------------------

;use run_adas405.pro to get the equilibrium ionisation balance
run_adas405, uid=uid, year=year,defyear=96,elem=elem, te=tein, dens=densin, frac=frac
nion=n_elements(frac.stage)
meta=fltarr(nion)
meta[*]=frac.ion[0,*]

;use run_adas406.pro to get the transient ionisation balance and
;run_adas405.pro to get the equilibrium ionisation balance for
;the same sets of Te,Ne (different from initial conditions
for t=0,ntime-1 do begin
for j=0,ndens-1 do begin
for i=0,nte-1 do begin
run_adas406, uid=uid, year=year, elem=elem, te=te[i], dens=dens[j], $
frac=frac, tint =time[t],meta=meta,defyear=96

tion[i]=frac.ion[0,is]
trec[i]=frac.ion[0,ir]
run_adas405, uid=uid, year=year, elem=elem, te=te[i], dens=dens[j], $
frac=fraceq,defyear=96

tioneq[i]=fraceq.ion[0,is]
treceq[i]=fraceq.ion[0,ir]
endfor
tdion[*,j]=tion[*]
tdrec[*,j]=trec[*]
ioneq[*,j]=tioneq[*]
receq[*,j]=treceq[*]
endfor
ion[*,*,t]=tdion[*,*]
rec[*,*,t]=tdrec[*,*]
endfor

;choose an electron density value
if not keyword_set(density) then begin
print,''
print,' Index Electron density values (cm-3):'
for d=0,ndens-1 do print,d,dens[d]
read,ii,prompt='Insert index for density (e.g. 3): '
ii=fix(ii)
endif

;---------------------------------------------------------
if keyword_set(norec) or keyword_set(exc) then begin
;---------------------------------------------------------
filenm=elem+strtrim(string(is),2)+'_'+strmid(strtrim(string(wlngth),2),0,6)+'_norec'
gnteq=ioneq*datas
for t=0,ntime-1 do gnt[*,*,t]=ion[*,*,t]*datas[*,*]
endif else begin
filenm=elem+strtrim(string(is),2)+'_'+strmid(strtrim(string(wlngth),2),0,6)
gnteq=ioneq*datas+receq*datar
for t=0,ntime-1 do gnt[*,*,t]=ion[*,*,t]*datas[*,*]+rec[*,*,t]*datar[*,*]
endelse

densi=dens[ii]
hh=reform(gnt[*,ii,*])
jj=fltarr(nte,ntime)
for g=0,ntime-1 do jj[*,g]=gnteq[*,ii]



;plot output
set_plot,'ps'
device, /isolatin1, font_index=8
device, bits=8, filename=filenm+'_'+strtrim(string(ii),2)+'_plot.ps', $
font_size = 14, xsize=16.0, ysize=18.0, $
yoffset=7.0, /color
device, /helvetica

!p.multi=[0,1,2]
surface,alog10(jj),alog10(tek),time,$
charsize=1.7,xtit='log(T!de!n) (K)',$
ytit='Time (s)',ztit='log(Contribution function)', $
tit='EQUILIBRIUM'


surface,alog10(hh),alog10(tek),time,$
charsize=1.7,xtit='log(T!de!n) (K)',$
ytit='Time (s)',ztit='log(Contribution function)', $
tit='TRANSIENT'

!p.multi=0

int=fix(nte/2)


plot,alog10(tek),alog10(hh[*,0]),xtit='log(Temperature) (K)',ytit='log(Contribution function)',line=1
for i=1,ntime-1 do oplot,alog10(tek),alog10(hh[*,i]),line=1
oplot,alog10(tek),alog10(gnteq[*,ii]),thick=3
xyouts,alog10(tek[int]),alog10(hh[int,0]),'t='+strmid(strtrim(string(min(time)),2),0,4)
xyouts,alog10(tek[int]),alog10(hh[int,ntime-1]),'t='+strmid(strtrim(string(max(time)),2),0,4)

device,/close
set_plot,'x'
!p.font=-1

print,''
print,' PLOT OUTPUT: '+filenm+'_'+strtrim(string(ii),2)+'_plot.ps'


;DA FARE ANCORA UN SAVE E UN TXT OUTPUT
;.save output
save,tekin,densin,tek,dens,time,gnteq,gnt,filename=filenm+'.save'

print,''
print,' SAVE OUTPUT: '+filenm+'.save'

;text output
openw,lun,filenm+'_'+strtrim(string(ii),2)+'_text.txt',/get_lun
printf,lun,' --------------------------------------- '
printf,lun,' INITIAL PLASMA CONDITIONS FOR TRANSIENT '
printf,lun,' Electron electron temperature (K)= ',tekin
printf,lun,' Electron electron density (cm-3)= ',densin
printf,lun,' --------------------------------------- '
printf,lun,' Selected final electron density (cm-3)= ',densi
printf,lun,' --------------------------------------- '
printf,lun,' --------------------------------------- '
printf,lun,' EQUILIBRIUM CONTRIBUTION FUNCTION FOR '+elem+strtrim(string(is),2)+' - '+strmid(strtrim(string(wlngth),2),0,6)+' Angstrom'
printf,lun,' Te (K) G(Te)'
for i=0,nte-1 do printf,lun,tek[i],gnteq[i,ii],format='(e9.2,3x,e9.2)'
printf,lun,' --------------------------------------- '
printf,lun,' --------------------------------------- '
printf,lun,' TRANSIENT CONTRIBUTION FUNCTION FOR '+elem+strtrim(string(is),2)+' - '+strmid(strtrim(string(wlngth),2),0,6)+' Angstrom'
printf,lun,'Te (K) time (s)'
printf,lun,''
printf,lun,' ',time,format='(a9,3x,'+strtrim(string(ntime),2)+'e9.2)'
printf,lun,''
for i=0,nte-1 do begin
printf,lun,tek[i],gnt[i,ii,*],format='(e9.2,3x,'+strtrim(string(ntime),2)+'e9.2)'
endfor


close,lun
free_lun,lun

print,''
print,' TEXT OUTPUT: '+filenm+'_'+strtrim(string(ii),2)+'_text.txt'


;stop
end