dustem_read_gas.pro
3.27 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
FUNCTION dustem_read_gas,file,silent=silent,help=help
;+
; NAME:
; dustem_read_gas
; PURPOSE:
; reads the file GAS.dat and returns the corresponding grain structure
; CATEGORY:
; Dustem
; CALLING SEQUENCE:
; st=dustem_read_gas(file,[/silent][,/help])
; INPUTS:
; file = file name to be read
; OPTIONAL INPUT PARAMETERS:
; None
; OUTPUTS:
; st: output structure
; OPTIONAL OUTPUT PARAMETERS:
; None
; ACCEPTED KEY-WORDS:
; silent = If set, keeps quiet
; help = If set, print this help
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; None
; RESTRICTIONS:
; The dustem idl wrapper must be installed
; PROCEDURE:
; None
; EXAMPLES
;
; MODIFICATION HISTORY:
; Written by J.-Ph. Bernard, N. Flagey, D. Paradis Jan 2007
; see evolution details on the dustem cvs maintained at CESR
; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-
;=== Format for DUSTEM_WHICH=WEB3p8
;# DUSTEM : Gas quantities from Ysard, Juvela & Verstraete (2011)
;#
;# Gas temperature (K), hydrogen density, H2 density, CR rate, G0 (if>0 supersedes GRAIN.DAT)
;# number of charge type (electrons, H+, C+, ...)
;# ion density (cm-3), mass (amu), charge, polarizability (A^3)
;# >>>>> 1st line is electron <<<<<<<
;#
; 3.5550e+03 1e-1 -1e0 5e-17 1.0000e+00 3
;-1e0 5.4858e-04 -1e0 0.00
;0e0 1.00794e0 1e0 0.67
; -1.3000e-05 12.0107e0 1e0 1.54
IF keyword_set(help) THEN BEGIN
doc_library,'dustem_read_gas'
st=0.
goto,the_end
ENDIF
CASE !dustem_which OF
'WEB3p8': BEGIN
;==read and count comments
openr,unit,file,/get_lun
str=''
comments=['']
first_char='#'
Ncomments=0L
WHILE first_char EQ '#' DO BEGIN
readf,unit,str
first_char=strmid(str,0,1)
comments=[comments,str]
IF first_char EQ '#' THEN Ncomments=Ncomments+1
ENDWHILE
comments=comments[1:*]
IF Ncomments NE 0 THEN comments=comments[0:Ncomments-1]
close,unit
free_lun,unit
st={file:file, $
Tgas:0.D0,nH:0.D0,nh2:0.D0,CR_rate:0.D0,G0:0.D0,n_charges:0L, $
charges_props:ptr_new(), $
comments:comments}
one_charge_st={density:0.D0,mass:0.D0,charge:0.D0,polarizability:0.D0}
;== pass comments
openr,unit,file,/get_lun
str=''
FOR i=0L,Ncomments-1 DO BEGIN
readf,unit,str
ENDFOR
;== read data
readf,unit,str
strv=str_sep(strcompress(strtrim(str,2)),' ') & ii=0L
st.Tgas=double(strv[ii]) & ii=ii+1
st.nH=double(strv[ii]) & ii=ii+1
st.nH2=double(strv[ii]) & ii=ii+1
st.CR_rate=double(strv[ii]) & ii=ii+1
st.G0=double(strv[ii]) & ii=ii+1
st.n_charges=long(strv[ii])
charges_st=replicate(one_charge_st,st.n_charges)
FOR i=0L,st.n_charges-1 DO BEGIN
readf,unit,str
strv=str_sep(strcompress(strtrim(str,2)),' ') & ii=0L
charges_st[i].density=double(strv[ii]) & ii=ii+1
charges_st[i].mass=double(strv[ii]) & ii=ii+1
charges_st[i].charge=double(strv[ii]) & ii=ii+1
charges_st[i].polarizability=double(strv[ii])
ENDFOR
close,unit
free_lun,unit
st.charges_props=ptr_new(charges_st)
END
'VERSTRAETE':BEGIN
;GAS.DAT did not exist in that version. Kept empty
st=0
END
ELSE:BEGIN
;GAS.DAT did not exist in that version. Kept empty
st=0
END
ENDCASE
;stop
the_end:
RETURN,st
END