OV-GSO-DC / dustem-wrapper_idl

Blame view

src/idl/dustem_read_grain.pro 5.62 KB
Edit Raw Normal View History
427f1205   Jean-Michel Glorian   version 4.2 merged 
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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
 
FUNCTION dustem_read_grain,file,silent=silent;,key_str=key_str,G0=G0,NH=NH

;=== Format for DUSTEM_WHICH=VERSTRAETE
;; # DUSTEM: definition of grain populations
;; # for each grain TYPE make sure you have the following files
;; # Q_TYPE.DAT in code_f90/les_QABS
;; # C_TYPE.DAT in code_f90/les_CAPA
;; # if option(MIX) MIX_TYPE.DAT in code_f90/les_DAT
;; # if option(POL) POL_TYPE.DAT in code_f90/les_DAT
;; # if keyword(SIZE) SIZE_TYPE.DAT in code_f90/les_DAT
;; # 
;; # run keywords
;; # G0 scaling factor for radiation field
;; # nr of grain types
;; # grain type -  Mdust/MH - rho(g/cm3) - amin(cm) - amax(cm) - alpha - nsize
;; #
;; SIZE
;; 1.00D+00
;; 3
;; PAH1    4.50D-04  2.25D+00  3.50D-08  7.22D-08  -3.50D+00   10  NONE
;; Gra     2.30D-03  2.25D+00  1.00D-07  5.00D-07  -3.50D+00   20  NONE
;; aSil    7.00D-03  3.30D+00  5.00D-07  1.10D-05  -3.50D+00   20  NONE
;; #

;=== Format for DUSTEM_WHICH=WEB3p8
;; # DUSTEM: definition of grain populations
;; # 
;; # run keywords 
;; # G0 scaling factor for radiation field
;; # grain type, nsize, type keywords, Mdust/MH, rho, amin, amax, alpha/a0 [, at, ac, gamma (ED)] [, au, zeta, eta (CV)]
;; # cgs units
;; #
;; sdist
;; 1.000000
;; PAH0              10 mix-logn           7.8000E-04  2.2400E+00  3.5000E-08  1.2000E-07  6.4000E-08  1.0000E-01
;; PAH1              10 mix-logn           7.8000E-04  2.2400E+00  3.5000E-08  1.2000E-07  6.4000E-08  1.0000E-01
;; amCBEx            15 logn         	   1.6500E-04  1.8100E+00  6.0000E-08  2.0000E-06  2.0000E-07  3.5000E-01
;; amCBEx  	     25 plaw-ed		   1.4500E-03  1.8100E+00  4.0000E-07  2.0000E-04 -2.8000E+00  1.5000E-05  1.5000E-05  2.0000E+00
;; aSil              25 plaw-ed   	   7.8000E-03  3.5000E+00  4.0000E-07  2.0000E-04 -3.4000E+00  2.0000E-05  2.0000E-05  2.0000E+00

;VG
defsysv, '!run_pol', 0.         ; Global flag to know if polarisation is being run
;DP
defsysv,'!run_tls',0.

;CASE getenv('DUSTEM_WHICH') OF
CASE !dustem_which OF
  'WEB3p8': BEGIN
    frmt='(A,D,D,D,D,D,I,A)'
    ;count # of lines
    openr,unit,file,/get_lun
    Nlines=0L
    str=''
    WHILE not eof(unit) DO BEGIN
      readf,unit,str
      Nlines=Nlines+1
    ENDWHILE
    close,unit
free_lun,unit
    ;== now read the file
    openr,unit,file,/get_lun
    ncurrent=0L
    REPEAT BEGIN
      readf,unit,str
      ncurrent=ncurrent+1
      first_char=strmid(str,0,1)
    ENDREP UNTIL first_char NE '#'
    key_str=str
;VG: we read G0, and NH if written beside G0
    ;readf,unit,G0,NH 
    readf,unit,str,format='(A100)'
    xx=strsplit(str,' ',/regex,/extr)
    G0 = float(xx(0))
    if n_elements(xx) gt 1 then NH=float(xx(1)) else NH=1.0
    ncurrent=ncurrent+1
    ;===
;    readf,unit,nspecies
    nspecies=nlines-ncurrent
;    stop
;; # grain type, nsize, type keywords, Mdust/MH, rho, amin, amax, alpha/a0 [, at, ac, gamma (ED)] [, au, zeta, eta (CV)]
    one_st={grain_type:'',nsize:0L,type_keywords:'',Mdust_o_MH:0.D0,rho:0.D0,amin:0.D0,amax:0.D0,alpha_o_a0:0.D0, $
            at:0.D0,ac:0.D0,gamma:0.D0, $
            au:0.D0,zeta:0.D0,eta:0.D0}
    st=replicate(one_st,nspecies)
    FOR i=0L,nspecies-1 DO BEGIN
      readf,unit,str
;VG : remove leading blanks
      ;strv=str_sep(strcompress(str),' ')
      strv=str_sep(strcompress(strtrim(str,2)),' ')
;VG
      Nstrv=n_elements(strv)
      ii=0L
      st(i).grain_type=strv(ii) & ii=ii+1
      st(i).nsize=long(strv(ii)) & ii=ii+1
      st(i).type_keywords=strv(ii) & ii=ii+1
;VG
      IF stregex(st(i).type_keywords, 'pol', /bool) THEN !run_pol = 1
;VG
;DP
      IF stregex(st(i).type_keywords, 'dtls', /bool) THEN !run_tls = 1

      st(i).Mdust_o_MH=double(strv(ii))/NH & ii=ii+1
      st(i).rho=double(strv(ii))  & ii=ii+1
      st(i).amin=double(strv(ii)) & ii=ii+1
      st(i).amax=double(strv(ii)) & ii=ii+1
      st(i).alpha_o_a0=double(strv(ii)) & ii=ii+1
      IF ii LE Nstrv-1 THEN BEGIN
        st(i).at=double(strv(ii)) & ii=ii+1
      ENDIF
      IF ii LE Nstrv-1 THEN BEGIN
        st(i).ac=double(strv(ii)) & ii=ii+1
      ENDIF
      IF ii LE Nstrv-1 THEN BEGIN
        st(i).gamma=double(strv(ii)) & ii=ii+1
      ENDIF
      IF ii LE Nstrv-1 THEN BEGIN
        st(i).au=double(strv(ii)) & ii=ii+1
      ENDIF
      IF ii LE Nstrv-1 THEN BEGIN
        st(i).zeta=double(strv(ii)) & ii=ii+1
      ENDIF
      IF ii LE Nstrv-1 THEN BEGIN
        st(i).eta=double(strv(ii)) & ii=ii+1
      ENDIF
    ENDFOR
    full_st={keywords:key_str,G0:G0,NH:NH,Ngrains:fix(nspecies),grains:st}
    ;full_st=st
    close,unit
    free_lun,unit
  END
  'VERSTRAETE':BEGIN
    frmt='(A,D,D,D,D,D,I,A)'
    openr,unit,file,/get_lun
    str=''
    REPEAT BEGIN
      readf,unit,str
      first_char=strmid(str,0,1)
    ENDREP UNTIL first_char NE '#'
    key_str=str
    readf,unit,G0
    readf,unit,nspecies
    one_st={type:'',propmass:0.D0,densgr:0.D0,tmin:0.D0,tmax:0.D0,alpha:0.D0,ndiscr:0,flag:''}
    st=replicate(one_st,nspecies)
    FOR i=0L,nspecies-1 DO BEGIN
      readf,unit,str
      strv=str_sep(strcompress(str),' ')
      st(i).type=strv(0)
      st(i).propmass=double(strv(1))
      st(i).densgr=double(strv(2))
      st(i).tmin=double(strv(3))
      st(i).tmax=double(strv(4))
      st(i).alpha=double(strv(5))
      st(i).ndiscr=fix(strv(6))
      st(i).flag=strv(7)
    ENDFOR
;    full_st={keywords:key_str,G0:G0,Ngrains:fix(nspecies),grains:st}
    full_st=st
    close,unit
    free_lun,unit
  END

  ELSE:BEGIN
    readcol,file,type,albmax,densgr,ngem,hydro,ioni,silent=silent
    nlines=n_elements(type)

    one_st={type:0,albmax:0.,densgr:0.,ngem:0,hydro:0.,ioni:0.}
    st=replicate(one_st,Nlines)
    st.type=type
    st.albmax=albmax
    st.densgr=densgr
    st.ngem=ngem
    st.hydro=hydro
    st.ioni=ioni
    full_st=st
  END
ENDCASE

RETURN,full_st


END