MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_thermal_emission.t
Go to the documentation of this file.
1! module mod_thermal_emission -- synthesize emission flux of some
2! thermal lines
3! EUV lines database:
4! 'He_II_304' 'Fe_IX_171' 'Fe_XXIV_193' 'Fe_XIV_211' 'Fe_XVI_335'
5! 'Fe_XVIII_94' 'Fe_XXI_131'
6! subroutines:
7! get_EUV: get local EUV emission intensity (for 1d, 2d and 3d)
8! get_SXR: get local Soft X-ray emission intensity (for 1d, 2d and 3d)
9
12 use mod_geometry
13 use mod_physics
14 use mod_comm_lib, only: mpistop
15
16 implicit none
17
18 integer :: n_aia
19 double precision :: t_aia(1:101)
20 double precision :: f_94(1:101),f_131(1:101),f_171(1:101)
21 double precision :: f_193(1:101),f_211(1:101),f_304(1:101)
22 double precision :: f_335(1:101)
23 integer :: n_iris
24 double precision :: t_iris(1:41)
25 double precision :: f_1354(1:41)
26 integer :: n_eis
27 double precision :: t_eis1(1:60),t_eis2(1:60)
28 double precision :: f_263(1:60),f_264(1:60),f_192(1:60),f_255(1:60)
29
30
31 double precision :: vec_xi1(1:3),vec_xi2(1:3),vec_los(1:3)
32
33 data n_aia / 101 /
34
35 data t_aia / 4. , 4.05, 4.1, 4.15, 4.2, 4.25, 4.3, 4.35, &
36 4.4, 4.45, 4.5, 4.55, 4.6, 4.65, 4.7, 4.75, &
37 4.8, 4.85, 4.9, 4.95, 5. , 5.05, 5.1, 5.15, &
38 5.2, 5.25, 5.3, 5.35, 5.4, 5.45, 5.5, 5.55, &
39 5.6, 5.65, 5.7, 5.75, 5.8, 5.85, 5.9, 5.95, &
40 6. , 6.05, 6.1, 6.15, 6.2, 6.25, 6.3, 6.35, &
41 6.4, 6.45, 6.5, 6.55, 6.6, 6.65, 6.7, 6.75, &
42 6.8, 6.85, 6.9, 6.95, 7. , 7.05, 7.1, 7.15, &
43 7.2, 7.25, 7.3, 7.35, 7.4, 7.45, 7.5, 7.55, &
44 7.6, 7.65, 7.7, 7.75, 7.8, 7.85, 7.9, 7.95, &
45 8. , 8.05, 8.1, 8.15, 8.2, 8.25, 8.3, 8.35, &
46 8.4, 8.45, 8.5, 8.55, 8.6, 8.65, 8.7, 8.75, &
47 8.8, 8.85, 8.9, 8.95, 9. /
48
49 data f_94 / 4.25022959d-37, 4.35880298d-36, 3.57054296d-35, 2.18175426d-34, &
50 8.97592571d-34, 2.68512961d-33, 7.49559346d-33, 2.11603751d-32, &
51 5.39752853d-32, 1.02935904d-31, 1.33822307d-31, 1.40884290d-31, &
52 1.54933156d-31, 2.07543102d-31, 3.42026227d-31, 6.31171444d-31, &
53 1.16559416d-30, 1.95360497d-30, 2.77818735d-30, 3.43552578d-30, &
54 4.04061803d-30, 4.75470982d-30, 5.65553769d-30, 6.70595782d-30, &
55 7.80680354d-30, 8.93247715d-30, 1.02618156d-29, 1.25979030d-29, &
56 1.88526483d-29, 3.62448572d-29, 7.50553279d-29, 1.42337571d-28, &
57 2.37912813d-28, 3.55232305d-28, 4.84985757d-28, 6.20662827d-28, &
58 7.66193687d-28, 9.30403645d-28, 1.10519802d-27, 1.25786927d-27, &
59 1.34362634d-27, 1.33185242d-27, 1.22302081d-27, 1.05677973d-27, &
60 9.23064720d-28, 8.78570994d-28, 8.02397416d-28, 5.87681142d-28, &
61 3.82272695d-28, 3.11492649d-28, 3.85736090d-28, 5.98893519d-28, &
62 9.57553548d-28, 1.46650267d-27, 2.10365847d-27, 2.79406671d-27, &
63 3.39420087d-27, 3.71077520d-27, 3.57296767d-27, 2.95114380d-27, &
64 2.02913103d-27, 1.13361825d-27, 5.13405629d-28, 2.01305089d-28, &
65 8.15781482d-29, 4.28366817d-29, 3.08701543d-29, 2.68693906d-29, &
66 2.51764203d-29, 2.41773103d-29, 2.33996083d-29, 2.26997246d-29, &
67 2.20316143d-29, 2.13810001d-29, 2.07424438d-29, 2.01149189d-29, &
68 1.94980213d-29, 1.88917920d-29, 1.82963583d-29, 1.77116920d-29, &
69 1.71374392d-29, 1.65740593d-29, 1.60214447d-29, 1.54803205d-29, &
70 1.49510777d-29, 1.44346818d-29, 1.39322305d-29, 1.34441897d-29, &
71 1.29713709d-29, 1.25132618d-29, 1.20686068d-29, 1.14226584d-29, &
72 1.09866413d-29, 1.05635524d-29, 1.01532444d-29, 9.75577134d-30, &
73 9.37102736d-30, 8.99873335d-30, 8.63860172d-30, 8.29051944d-30, &
74 7.95414793d-30 /
75
76 data f_131 / 3.18403601d-37, 3.22254703d-36, 2.61657920d-35, &
77 1.59575286d-34, 6.65779556d-34, 2.07015132d-33, &
78 6.05768615d-33, 1.76074833d-32, 4.52633001d-32, &
79 8.57121883d-32, 1.09184271d-31, 1.10207963d-31, &
80 1.11371658d-31, 1.29105226d-31, 1.80385897d-31, &
81 3.27295431d-31, 8.92002136d-31, 3.15214579d-30, &
82 9.73440787d-30, 2.22709702d-29, 4.01788984d-29, &
83 6.27471832d-29, 8.91764995d-29, 1.18725647d-28, &
84 1.52888040d-28, 2.05082946d-28, 3.47651873d-28, &
85 8.80482184d-28, 2.66533063d-27, 7.05805149d-27, &
86 1.46072515d-26, 2.45282476d-26, 3.55303726d-26, &
87 4.59075911d-26, 5.36503515d-26, 5.68444094d-26, &
88 5.47222296d-26, 4.81119761d-26, 3.85959059d-26, &
89 2.80383406d-26, 1.83977650d-26, 1.11182849d-26, &
90 6.50748885d-27, 3.96843481d-27, 2.61876319d-27, &
91 1.85525324d-27, 1.39717024d-27, 1.11504283d-27, &
92 9.38169611d-28, 8.24801234d-28, 7.43331919d-28, &
93 6.74537063d-28, 6.14495760d-28, 5.70805277d-28, &
94 5.61219786d-28, 6.31981777d-28, 9.19747307d-28, &
95 1.76795732d-27, 3.77985446d-27, 7.43166191d-27, &
96 1.19785603d-26, 1.48234676d-26, 1.36673114d-26, &
97 9.61047146d-27, 5.61209353d-27, 3.04779780d-27, &
98 1.69378976d-27, 1.02113491d-27, 6.82223774d-28, &
99 5.02099099d-28, 3.99377760d-28, 3.36279037d-28, &
100 2.94767378d-28, 2.65740865d-28, 2.44396277d-28, &
101 2.28003967d-28, 2.14941419d-28, 2.04178995d-28, &
102 1.95031045d-28, 1.87011994d-28, 1.79777869d-28, &
103 1.73093957d-28, 1.66795789d-28, 1.60785455d-28, &
104 1.55002399d-28, 1.49418229d-28, 1.44022426d-28, &
105 1.38807103d-28, 1.33772767d-28, 1.28908404d-28, &
106 1.24196208d-28, 1.17437501d-28, 1.12854330d-28, &
107 1.08410498d-28, 1.04112003d-28, 9.99529904d-29, &
108 9.59358806d-29, 9.20512291d-29, 8.83009123d-29, &
109 8.46817043d-29, 8.11921928d-29 /
110
111 data f_171 / 2.98015581d-42, 1.24696230d-40, 3.37614652d-39, 5.64103034d-38, &
112 5.20550266d-37, 2.77785939d-36, 1.16283616d-35, 6.50007689d-35, &
113 9.96177399d-34, 1.89586076d-32, 2.10982799d-31, 1.36946479d-30, &
114 6.27396553d-30, 2.29955134d-29, 7.13430211d-29, 1.91024282d-28, &
115 4.35358848d-28, 7.94807808d-28, 1.07431875d-27, 1.08399488d-27, &
116 9.16212938d-28, 7.34715770d-28, 6.59246382d-28, 9.13541375d-28, &
117 2.05939035d-27, 5.08206555d-27, 1.10148083d-26, 2.01884662d-26, &
118 3.13578384d-26, 4.14367719d-26, 5.36067711d-26, 8.74170213d-26, &
119 1.64161233d-25, 2.94587860d-25, 4.76298332d-25, 6.91765639d-25, &
120 9.08825111d-25, 1.08496183d-24, 1.17440114d-24, 1.13943939d-24, &
121 9.71696981d-25, 7.09593688d-25, 4.31376399d-25, 2.12708486d-25, &
122 8.47429567d-26, 3.17608104d-26, 1.95898842d-26, 1.98064242d-26, &
123 1.67706555d-26, 8.99126003d-27, 3.29773878d-27, 1.28896127d-27, &
124 8.51169698d-28, 7.53520167d-28, 6.18268143d-28, 4.30034650d-28, &
125 2.78152409d-28, 1.95437088d-28, 1.65896278d-28, 1.68740181d-28, &
126 1.76054383d-28, 1.63978419d-28, 1.32880591d-28, 1.00833205d-28, &
127 7.82252806d-29, 6.36181741d-29, 5.34633869d-29, 4.58013864d-29, &
128 3.97833422d-29, 3.49414760d-29, 3.09790940d-29, 2.76786227d-29, &
129 2.48806269d-29, 2.24823367d-29, 2.04016653d-29, 1.85977413d-29, &
130 1.70367499d-29, 1.56966125d-29, 1.45570643d-29, 1.35964565d-29, &
131 1.27879263d-29, 1.21016980d-29, 1.15132499d-29, 1.09959628d-29, &
132 1.05307482d-29, 1.01040261d-29, 9.70657096d-30, 9.33214234d-30, &
133 8.97689427d-30, 8.63761192d-30, 8.31149879d-30, 7.85162401d-30, &
134 7.53828281d-30, 7.23559452d-30, 6.94341530d-30, 6.66137038d-30, &
135 6.38929156d-30, 6.12669083d-30, 5.87346434d-30, 5.62943622d-30, &
136 5.39435202d-30 /
137
138 data f_193 / 6.40066486d-32, 4.92737300d-31, 2.95342934d-30, 1.28061594d-29, &
139 3.47747667d-29, 5.88554792d-29, 7.72171179d-29, 9.75609282d-29, &
140 1.34318963d-28, 1.96252638d-28, 2.70163878d-28, 3.63192965d-28, &
141 5.28087341d-28, 8.37821446d-28, 1.39089159d-27, 2.31749718d-27, &
142 3.77510689d-27, 5.85198594d-27, 8.26021568d-27, 1.04870405d-26, &
143 1.25209374d-26, 1.47406787d-26, 1.77174067d-26, 2.24098537d-26, &
144 3.05926105d-26, 4.50018853d-26, 6.84720216d-26, 1.00595861d-25, &
145 1.30759222d-25, 1.36481773d-25, 1.15943558d-25, 1.01467304d-25, &
146 1.04092532d-25, 1.15071251d-25, 1.27416033d-25, 1.38463476d-25, &
147 1.47882726d-25, 1.57041238d-25, 1.69786224d-25, 1.94970397d-25, &
148 2.50332918d-25, 3.58321431d-25, 5.18061550d-25, 6.60405549d-25, &
149 6.64085365d-25, 4.83825816d-25, 2.40545020d-25, 8.59534098d-26, &
150 2.90920638d-26, 1.33204845d-26, 9.03933926d-27, 7.78910836d-27, &
151 7.29342321d-27, 7.40267022d-27, 8.05279981d-27, 8.13829291d-27, &
152 6.92634262d-27, 5.12521880d-27, 3.59527615d-27, 2.69617560d-27, &
153 2.84432713d-27, 5.06697306d-27, 1.01281903d-26, 1.63526978d-26, &
154 2.06759342d-26, 2.19482312d-26, 2.10050611d-26, 1.89837248d-26, &
155 1.66347131d-26, 1.43071097d-26, 1.21518419d-26, 1.02078343d-26, &
156 8.46936184d-27, 6.93015742d-27, 5.56973237d-27, 4.38951754d-27, &
157 3.38456457d-27, 2.55309556d-27, 1.88904224d-27, 1.38057546d-27, &
158 1.00718330d-27, 7.43581116d-28, 5.63562931d-28, 4.43359435d-28, &
159 3.63923535d-28, 3.11248143d-28, 2.75586846d-28, 2.50672237d-28, &
160 2.32419348d-28, 2.18325682d-28, 2.06834486d-28, 1.93497044d-28, &
161 1.84540751d-28, 1.76356504d-28, 1.68741425d-28, 1.61566157d-28, &
162 1.54754523d-28, 1.48249410d-28, 1.42020176d-28, 1.36045230d-28, &
163 1.30307965d-28 /
164
165 data f_211 / 4.74439912d-42, 1.95251522d-40, 5.19700194d-39, 8.53120166d-38, &
166 7.72745727d-37, 4.04158559d-36, 1.64853511d-35, 8.56295439d-35, &
167 1.17529722d-33, 2.16867729d-32, 2.40472264d-31, 1.56418133d-30, &
168 7.20032889d-30, 2.65838271d-29, 8.33196904d-29, 2.26128236d-28, &
169 5.24295811d-28, 9.77791121d-28, 1.35913489d-27, 1.43957785d-27, &
170 1.37591544d-27, 1.49029886d-27, 2.06183401d-27, 3.31440622d-27, &
171 5.42497318d-27, 8.41100374d-27, 1.17941366d-26, 1.49269794d-26, &
172 1.71506074d-26, 1.71266353d-26, 1.51434781d-26, 1.36766622d-26, &
173 1.33483562d-26, 1.36834518d-26, 1.45829002d-26, 1.62575306d-26, &
174 1.88773347d-26, 2.22026986d-26, 2.54930499d-26, 2.80758138d-26, &
175 3.06176409d-26, 3.62799792d-26, 5.13226109d-26, 8.46260744d-26, &
176 1.38486586d-25, 1.86192535d-25, 1.78007934d-25, 1.16548409d-25, &
177 5.89293257d-26, 2.69952884d-26, 1.24891081d-26, 6.41273176d-27, &
178 4.08282914d-27, 3.26463328d-27, 2.76230280d-27, 2.08986882d-27, &
179 1.37658470d-27, 8.48489381d-28, 5.19304217d-28, 3.19312514d-28, &
180 2.02968197d-28, 1.50171666d-28, 1.39164218d-28, 1.42448821d-28, &
181 1.41714519d-28, 1.33341059d-28, 1.20759270d-28, 1.07259692d-28, &
182 9.44895400d-29, 8.29030041d-29, 7.25440631d-29, 6.33479483d-29, &
183 5.51563757d-29, 4.79002469d-29, 4.14990482d-29, 3.59384972d-29, &
184 3.12010860d-29, 2.72624742d-29, 2.40734791d-29, 2.15543565d-29, &
185 1.95921688d-29, 1.80682882d-29, 1.68695662d-29, 1.59020936d-29, &
186 1.50940886d-29, 1.43956179d-29, 1.37731622d-29, 1.32049043d-29, &
187 1.26771875d-29, 1.21803879d-29, 1.17074716d-29, 1.10507836d-29, &
188 1.06022834d-29, 1.01703080d-29, 9.75436986d-30, 9.35349257d-30, &
189 8.96744546d-30, 8.59527489d-30, 8.23678940d-30, 7.89144480d-30, &
190 7.55891138d-30 /
191
192 data f_304 / 3.62695850d-32, 2.79969087d-31, 1.68340584d-30, 7.32681440d-30, &
193 1.99967770d-29, 3.41296785d-29, 4.55409104d-29, 5.94994635d-29, &
194 8.59864963d-29, 1.39787633d-28, 3.17701965d-28, 1.14474920d-27, &
195 4.44845958d-27, 1.54785841d-26, 4.70265345d-26, 1.24524365d-25, &
196 2.81535352d-25, 5.10093666d-25, 6.83545307d-25, 6.82110329d-25, &
197 5.66886188d-25, 4.36205513d-25, 3.29265688d-25, 2.49802368d-25, &
198 1.92527113d-25, 1.51058572d-25, 1.20596047d-25, 9.76884267d-26, &
199 7.89979266d-26, 6.18224289d-26, 4.67298332d-26, 3.57934505d-26, &
200 2.84535785d-26, 2.32853022d-26, 1.95228514d-26, 1.67880071d-26, &
201 1.47608785d-26, 1.32199691d-26, 1.20070960d-26, 1.09378177d-26, &
202 1.00031730d-26, 9.62434001d-27, 1.05063954d-26, 1.27267143d-26, &
203 1.45923057d-26, 1.36746707d-26, 1.03466970d-26, 6.97647829d-27, &
204 4.63141039d-27, 3.19031994d-27, 2.33373613d-27, 1.81589079d-27, &
205 1.48446917d-27, 1.26611478d-27, 1.12617468d-27, 1.03625148d-27, &
206 9.61400595d-28, 8.79016231d-28, 7.82612130d-28, 6.73762960d-28, &
207 5.59717956d-28, 4.53010243d-28, 3.65712196d-28, 3.00958686d-28, &
208 2.54011502d-28, 2.18102277d-28, 1.88736437d-28, 1.63817539d-28, &
209 1.42283147d-28, 1.23631916d-28, 1.07526003d-28, 9.36797928d-29, &
210 8.18565660d-29, 7.18152734d-29, 6.32523238d-29, 5.59513985d-29, &
211 4.96614048d-29, 4.42518826d-29, 3.95487628d-29, 3.54690294d-29, &
212 3.18953930d-29, 2.87720933d-29, 2.60186750d-29, 2.36011522d-29, &
213 2.14717806d-29, 1.95905217d-29, 1.79287981d-29, 1.64562262d-29, &
214 1.51489425d-29, 1.39876064d-29, 1.29496850d-29, 1.18665438d-29, &
215 1.10240474d-29, 1.02643099d-29, 9.57780996d-30, 8.95465151d-30, &
216 8.38950190d-30, 7.87283711d-30, 7.40136507d-30, 6.96804279d-30, &
217 6.56945323d-30 /
218
219 data f_335 / 2.46882661d-32, 1.89476632d-31, 1.13216502d-30, 4.89532008d-30, &
220 1.32745970d-29, 2.25390335d-29, 3.00511672d-29, 3.96035934d-29, &
221 5.77977656d-29, 8.58600736d-29, 1.14083000d-28, 1.48644411d-28, &
222 2.15788823d-28, 3.51628877d-28, 6.12200698d-28, 1.08184987d-27, &
223 1.85590697d-27, 2.91679107d-27, 3.94405396d-27, 4.63610680d-27, &
224 5.13824456d-27, 5.66602209d-27, 6.30009232d-27, 7.03422868d-27, &
225 7.77973918d-27, 8.32371831d-27, 8.56724316d-27, 8.62601374d-27, &
226 8.13308844d-27, 6.53188216d-27, 4.55197029d-27, 3.57590087d-27, &
227 3.59571707d-27, 4.03502770d-27, 4.54366411d-27, 4.96914990d-27, &
228 5.24601170d-27, 5.39979250d-27, 5.43023669d-27, 5.26235042d-27, &
229 4.91585495d-27, 4.52628362d-27, 4.13385020d-27, 3.67538967d-27, &
230 3.39939742d-27, 3.81284533d-27, 5.02332701d-27, 6.19438602d-27, &
231 6.49613071d-27, 6.04010475d-27, 5.24664275d-27, 4.37225997d-27, &
232 3.52957182d-27, 2.76212276d-27, 2.08473158d-27, 1.50850518d-27, &
233 1.04602472d-27, 7.13091243d-28, 5.34289645d-28, 5.21079581d-28, &
234 6.22246365d-28, 6.99555864d-28, 6.29665489d-28, 4.45077026d-28, &
235 2.67046793d-28, 1.52774686d-28, 9.18061770d-29, 6.09116074d-29, &
236 4.48562572d-29, 3.59463696d-29, 3.05820218d-29, 2.70766652d-29, &
237 2.46144034d-29, 2.27758450d-29, 2.13331183d-29, 2.01537836d-29, &
238 1.91566180d-29, 1.82893912d-29, 1.75167748d-29, 1.68136168d-29, &
239 1.61615595d-29, 1.55481846d-29, 1.49643236d-29, 1.44046656d-29, &
240 1.38657085d-29, 1.33459068d-29, 1.28447380d-29, 1.23615682d-29, &
241 1.18963296d-29, 1.14478976d-29, 1.10146637d-29, 1.04039479d-29, &
242 9.98611410d-30, 9.58205147d-30, 9.19202009d-30, 8.81551313d-30, &
243 8.45252127d-30, 8.10224764d-30, 7.76469090d-30, 7.43954323d-30, &
244 7.12653873d-30 /
245
246
247 data n_iris / 41 /
248
249 data t_iris / 4. , 4.1 , 4.2 , 4.3 , 4.40000001, &
250 4.50000001, 4.60000001, 4.70000001, 4.80000001, 4.90000001, &
251 5.00000001, 5.10000002, 5.20000002, 5.30000002, 5.40000002, &
252 5.50000002, 5.60000002, 5.70000003, 5.80000003, 5.90000003, &
253 6.00000003, 6.10000003, 6.20000003, 6.30000003, 6.40000004, &
254 6.50000004, 6.60000004, 6.70000004, 6.80000004, 6.90000004, &
255 7.00000004, 7.10000005, 7.20000005, 7.30000005, 7.40000005, &
256 7.50000005, 7.60000005, 7.70000006, 7.80000006, 7.90000006, &
257 8.00000006 /
258
259 data f_1354 / 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
260 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
261 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
262 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
263 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
264 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 1.09503647d-39, &
265 5.47214550d-36, 2.42433983d-33, 2.75295034d-31, 1.21929718d-29, &
266 2.48392125d-28, 2.33268145d-27, 8.68623633d-27, 1.00166284d-26, &
267 3.63126633d-27, 7.45174807d-28, 1.38224064d-28, 2.69270994d-29, &
268 5.53314977d-30, 1.15313092d-30, 2.34195788d-31, 4.48242942d-32, &
269 7.94976380d-33 /
270
271
272 data n_eis / 60 /
273
274 data t_eis1 / 1.99526231d+05, 2.23872114d+05, 2.51188643d+05, 2.81838293d+05, &
275 3.16227766d+05, 3.54813389d+05, 3.98107171d+05, 4.46683592d+05, &
276 5.01187234d+05, 5.62341325d+05, 6.30957344d+05, 7.07945784d+05, &
277 7.94328235d+05, 8.91250938d+05, 1.00000000d+06, 1.12201845d+06, &
278 1.25892541d+06, 1.41253754d+06, 1.58489319d+06, 1.77827941d+06, &
279 1.99526231d+06, 2.23872114d+06, 2.51188643d+06, 2.81838293d+06, &
280 3.16227766d+06, 3.54813389d+06, 3.98107171d+06, 4.46683592d+06, &
281 5.01187234d+06, 5.62341325d+06, 6.30957344d+06, 7.07945784d+06, &
282 7.94328235d+06, 8.91250938d+06, 1.00000000d+07, 1.12201845d+07, &
283 1.25892541d+07, 1.41253754d+07, 1.58489319d+07, 1.77827941d+07, &
284 1.99526231d+07, 2.23872114d+07, 2.51188643d+07, 2.81838293d+07, &
285 3.16227766d+07, 3.54813389d+07, 3.98107171d+07, 4.46683592d+07, &
286 5.01187234d+07, 5.62341325d+07, 6.30957344d+07, 7.07945784d+07, &
287 7.94328235d+07, 8.91250938d+07, 1.00000000d+08, 1.12201845d+08, &
288 1.25892541d+08, 1.41253754d+08, 1.58489319d+08, 1.77827941d+08 /
289
290 data t_eis2 / 1.99526231d+06, 2.23872114d+06, 2.51188643d+06, 2.81838293d+06, &
291 3.16227766d+06, 3.54813389d+06, 3.98107171d+06, 4.46683592d+06, &
292 5.01187234d+06, 5.62341325d+06, 6.30957344d+06, 7.07945784d+06, &
293 7.94328235d+06, 8.91250938d+06, 1.00000000d+07, 1.12201845d+07, &
294 1.25892541d+07, 1.41253754d+07, 1.58489319d+07, 1.77827941d+07, &
295 1.99526231d+07, 2.23872114d+07, 2.51188643d+07, 2.81838293d+07, &
296 3.16227766d+07, 3.54813389d+07, 3.98107171d+07, 4.46683592d+07, &
297 5.01187234d+07, 5.62341325d+07, 6.30957344d+07, 7.07945784d+07, &
298 7.94328235d+07, 8.91250938d+07, 1.00000000d+08, 1.12201845d+08, &
299 1.25892541d+08, 1.41253754d+08, 1.58489319d+08, 1.77827941d+08, &
300 1.99526231d+08, 2.23872114d+08, 2.51188643d+08, 2.81838293d+08, &
301 3.16227766d+08, 3.54813389d+08, 3.98107171d+08, 4.46683592d+08, &
302 5.01187234d+08, 5.62341325d+08, 6.30957344d+08, 7.07945784d+08, &
303 7.94328235d+08, 8.91250938d+08, 1.00000000d+09, 1.12201845d+09, &
304 1.25892541d+09, 1.41253754d+09, 1.58489319d+09, 1.77827941d+09 /
305
306 data f_263 / 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
307 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
308 0.00000000d+00, 4.46454917d-45, 3.26774829d-42, 1.25292566d-39, &
309 2.66922338d-37, 3.28497742d-35, 2.38677554d-33, 1.03937729d-31, &
310 2.75075687d-30, 4.47961733d-29, 4.46729177d-28, 2.64862689d-27, &
311 8.90863800d-27, 1.72437548d-26, 2.22217752d-26, 2.27999477d-26, &
312 2.08264363d-26, 1.78226687d-26, 1.45821699d-26, 1.14675379d-26, &
313 8.63082492d-27, 6.15925429d-27, 4.11252514d-27, 2.51530564d-27, &
314 1.37090986d-27, 6.42443134d-28, 2.48392636d-28, 7.59187874d-29, &
315 1.77852938d-29, 3.23945221d-30, 4.90533903d-31, 6.75458158d-32, &
316 9.06878868d-33, 1.23927474d-33, 1.75769395d-34, 2.60710914d-35, &
317 4.04318030d-36, 6.53500581d-37, 1.09365022d-37, 1.88383322d-38, &
318 3.31425233d-39, 5.90964084d-40, 1.06147549d-40, 1.90706170d-41, &
319 3.41331584d-42, 6.07310718d-43, 1.07364738d-43, 1.89085498d-44, &
320 3.32598922d-45, 5.87125640d-46, 0.00000000d+00, 0.00000000d+00 /
321
322 data f_264 / 0.00000000d+00, 2.81670057d-46, 1.28007268d-43, 2.54586603d-41, &
323 2.67887256d-39, 1.68413285d-37, 6.85702304d-36, 1.91797284d-34, &
324 3.84675839d-33, 5.69939170d-32, 6.36224608d-31, 5.39176489d-30, &
325 3.45478458d-29, 1.64848693d-28, 5.71476364d-28, 1.39909997d-27, &
326 2.37743056d-27, 2.86712530d-27, 2.65206348d-27, 2.07175767d-27, &
327 1.47866767d-27, 1.01087374d-27, 6.79605811d-28, 4.54746770d-28, &
328 3.04351751d-28, 2.03639149d-28, 1.35940991d-28, 9.01451939d-29, &
329 5.91289972d-29, 3.81821178d-29, 2.41434696d-29, 1.48871220d-29, &
330 8.93362094d-30, 5.21097445d-30, 2.95964719d-30, 1.64278748d-30, &
331 8.95571660d-31, 4.82096011d-31, 2.57390991d-31, 1.36821781d-31, &
332 7.27136350d-32, 3.87019426d-32, 2.06883430d-32, 1.11228884d-32, &
333 6.01883313d-33, 3.27790676d-33, 1.79805012d-33, 9.93085346d-34, &
334 5.52139556d-34, 3.08881387d-34, 1.73890315d-34, 9.84434964d-35, &
335 5.60603378d-35, 3.20626492d-35, 1.84111068d-35, 0.00000000d+00, &
336 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00 /
337
338 data f_192 / 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 4.35772105d-44, &
339 1.26162319d-41, 1.97471205d-39, 1.83409019d-37, 1.08206288d-35, &
340 4.27914363d-34, 1.17943846d-32, 2.32565755d-31, 3.33087991d-30, &
341 3.47013260d-29, 2.60375866d-28, 1.37737127d-27, 5.01053913d-27, &
342 1.23479810d-26, 2.11310542d-26, 2.71831513d-26, 2.89851163d-26, &
343 2.77312376d-26, 2.50025229d-26, 2.18323661d-26, 1.86980322d-26, &
344 1.58035034d-26, 1.31985651d-26, 1.08733133d-26, 8.81804906d-27, &
345 7.00417973d-27, 5.43356567d-27, 4.09857884d-27, 2.99651764d-27, &
346 2.11902962d-27, 1.45014127d-27, 9.62291023d-28, 6.21548647d-28, &
347 3.92807578d-28, 2.44230375d-28, 1.50167782d-28, 9.17611405d-29, &
348 5.58707641d-29, 3.40570915d-29, 2.08030862d-29, 1.27588676d-29, &
349 7.86535588d-30, 4.87646151d-30, 3.03888897d-30, 1.90578649d-30, &
350 1.20195947d-30, 7.61955060d-31, 4.85602199d-31, 3.11049969d-31, &
351 2.00087065d-31, 1.29223740d-31, 8.37422008d-32, 0.00000000d+00, &
352 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00 /
353
354 data f_255 / 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 1.76014287d-44, &
355 5.07057938d-42, 7.90473970d-40, 7.31852999d-38, 4.30709255d-36, &
356 1.70009061d-34, 4.67925160d-33, 9.21703546d-32, 1.31918676d-30, &
357 1.37393161d-29, 1.03102379d-28, 5.45694018d-28, 1.98699648d-27, &
358 4.90346776d-27, 8.40524725d-27, 1.08321456d-26, 1.15714525d-26, &
359 1.10905152d-26, 1.00155023d-26, 8.75799161d-27, 7.50935839d-27, &
360 6.35253533d-27, 5.30919268d-27, 4.37669455d-27, 3.55185164d-27, &
361 2.82347055d-27, 2.19257595d-27, 1.65589541d-27, 1.21224987d-27, &
362 8.58395132d-28, 5.88163935d-28, 3.90721447d-28, 2.52593407d-28, &
363 1.59739995d-28, 9.93802874d-29, 6.11343388d-29, 3.73711135d-29, &
364 2.27618743d-29, 1.38793199d-29, 8.48060787d-30, 5.20305940d-30, &
365 3.20867365d-30, 1.99011277d-30, 1.24064551d-30, 7.78310544d-31, &
366 4.91013681d-31, 3.11338381d-31, 1.98451675d-31, 1.27135460d-31, &
367 8.17917486d-32, 5.28280497d-32, 3.42357159d-32, 0.00000000d+00, &
368 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00 /
369
370 abstract interface
371 subroutine get_subr1(w,x,ixI^L,ixO^L,res)
373 integer, intent(in) :: ixI^L, ixO^L
374 double precision, intent(in) :: w(ixI^S,nw)
375 double precision, intent(in) :: x(ixI^S,1:ndim)
376 double precision, intent(out):: res(ixI^S)
377 end subroutine get_subr1
378
379 end interface
380
382
383 procedure(get_subr1), pointer, nopass :: get_rho => null()
384 procedure(get_subr1), pointer, nopass :: get_pthermal => null()
385 procedure(get_subr1), pointer, nopass :: get_var_rfactor => null()
386
387 end type te_fluid
388
389
390 contains
391
392 subroutine get_line_info(wl,ion,mass,logTe,line_center,spatial_px,spectral_px,sigma_PSF,width_slit)
393 ! get information of the spectral line
394 ! wl: wavelength
395 ! mass: ion mass, unit -- proton mass
396 ! logTe: peak temperature of emission line in logarithm
397 ! line_center: center wavelength of emission line, unit -- Angstrom (0.1 nm)
398 ! spatial_px: pixel size in space of instrument (for image), unit -- arcsec
399 ! spectral_px: pixel size in wagelength of instrument (for spectrum), unit -- Angstrom
400 ! sigma_PSF: width of point spread function core (for instrument), unit -- pixel
401 ! width_slit: width of slit for spectrograph, unit -- arcsec
403
404 integer, intent(in) :: wl
405 integer, intent(out) :: mass
406 character(len=30), intent(out) :: ion
407 double precision, intent(out) :: logTe,line_center,spatial_px,spectral_px
408 double precision, intent(out) :: sigma_PSF,width_slit
409
410 select case(wl)
411 case(304)
412 ion='He II'
413 mass=4
414 logte=4.7d0
415 line_center=303.8d0
416 spatial_px=0.6d0
417 spectral_px=0.02d0
418 sigma_psf=0.895d0
419 width_slit=0.6d0
420 case(171)
421 ion='Fe IX'
422 mass=56
423 logte=5.8d0
424 line_center=171.1d0
425 spatial_px=0.6d0
426 spectral_px=0.02d0
427 sigma_psf=1.019d0
428 width_slit=0.6d0
429 case(193)
430 ion='Fe XXIV'
431 mass=56
432 logte=7.3d0
433 line_center=193.5d0
434 spatial_px=0.6d0
435 spectral_px=0.02d0
436 sigma_psf=0.813d0
437 width_slit=0.6d0
438 case(211)
439 ion='Fe XIV'
440 mass=56
441 logte=6.3d0
442 line_center=211.3d0
443 spatial_px=0.6d0
444 spectral_px=0.02d0
445 sigma_psf=0.913d0
446 width_slit=0.6d0
447 case(335)
448 ion='Fe XVI'
449 mass=56
450 logte=6.4d0
451 line_center=335.4d0
452 spatial_px=0.6d0
453 spectral_px=0.02d0
454 sigma_psf=1.019d0
455 width_slit=0.6d0
456 case(94)
457 ion='Fe XVIII'
458 mass=56
459 logte=6.8d0
460 line_center=93.9d0
461 spatial_px=0.6d0
462 spectral_px=0.02d0
463 sigma_psf=1.025d0
464 width_slit=0.6d0
465 case(131)
466 ion='Fe XXI'
467 mass=56
468 logte=7.0d0
469 line_center=131.0d0
470 spatial_px=0.6d0
471 spectral_px=0.02d0
472 sigma_psf=0.984d0
473 width_slit=0.6d0
474 case(1354)
475 ion='Fe XXI'
476 mass=56
477 logte=7.0d0
478 line_center=1354.1d0
479 spatial_px=0.1663d0
480 spectral_px=12.98d-3
481 sigma_psf=1.d0
482 width_slit=0.33d0
483 case(263)
484 ion='Fe XVI'
485 mass=56
486 logte=6.4d0
487 line_center=262.976d0
488 spatial_px=1.d0
489 spectral_px=22.d-3
490 sigma_psf=1.d0
491 width_slit=2.d0
492 case(264)
493 ion='Fe XXIII'
494 mass=56
495 logte=7.1d0
496 line_center=263.765d0
497 spatial_px=1.d0
498 spectral_px=22.d-3
499 sigma_psf=1.d0
500 width_slit=2.d0
501 case(192)
502 ion='Fe XXIV'
503 mass=56
504 logte=7.2d0
505 line_center=192.028d0
506 spatial_px=1.d0
507 spectral_px=22.d-3
508 sigma_psf=1.d0
509 width_slit=2.d0
510 case(255)
511 ion='Fe XXIV'
512 mass=56
513 logte=7.2d0
514 line_center=255.113d0
515 spatial_px=1.d0
516 spectral_px=22.d-3
517 sigma_psf=1.d0
518 width_slit=2.d0
519 case default
520 call mpistop("No information about this line")
521 end select
522
523 spatial_px=spatial_px/instrument_resolution_factor
524 end subroutine get_line_info
525
526 subroutine get_euv(wl,ixI^L,ixO^L,w,x,fl,flux)
527 ! calculate the local emission intensity of given EUV line (optically thin)
528 ! wavelength is the wave length of the emission line
529 ! unit [DN cm^-1 s^-1 pixel^-1]
530 ! ingrate flux along line of sight: DN s^-1 pixel^-1
532
533 integer, intent(in) :: wl
534 integer, intent(in) :: ixI^L, ixO^L
535 double precision, intent(in) :: x(ixI^S,1:ndim)
536 double precision, intent(in) :: w(ixI^S,1:nw)
537 type(te_fluid), intent(in) :: fl
538 double precision, intent(out) :: flux(ixI^S)
539
540 integer :: n_table
541 double precision, allocatable :: t_table(:),f_table(:)
542 integer :: ix^D,iTt,i
543 double precision :: pth(ixI^S),Te(ixI^S),Ne(ixI^S)
544 double precision :: logT,logGT,GT
545
546 ! selecting emission table
547 select case(wl)
548 case(94)
549 n_table=n_aia
550 allocate(t_table(1:n_table))
551 allocate(f_table(1:n_table))
552 t_table(1:n_table)=t_aia(1:n_aia)
553 f_table(1:n_table)=f_94(1:n_aia)
554 case(131)
555 n_table=n_aia
556 allocate(t_table(1:n_table))
557 allocate(f_table(1:n_table))
558 t_table(1:n_table)=t_aia(1:n_aia)
559 f_table(1:n_table)=f_131(1:n_aia)
560 case(171)
561 n_table=n_aia
562 allocate(t_table(1:n_table))
563 allocate(f_table(1:n_table))
564 t_table(1:n_table)=t_aia(1:n_aia)
565 f_table(1:n_table)=f_171(1:n_aia)
566 case(193)
567 n_table=n_aia
568 allocate(t_table(1:n_table))
569 allocate(f_table(1:n_table))
570 t_table(1:n_table)=t_aia(1:n_aia)
571 f_table(1:n_table)=f_193(1:n_aia)
572 case(211)
573 n_table=n_aia
574 allocate(t_table(1:n_table))
575 allocate(f_table(1:n_table))
576 t_table(1:n_table)=t_aia(1:n_aia)
577 f_table(1:n_table)=f_211(1:n_aia)
578 case(304)
579 n_table=n_aia
580 allocate(t_table(1:n_table))
581 allocate(f_table(1:n_table))
582 t_table(1:n_table)=t_aia(1:n_aia)
583 f_table(1:n_table)=f_304(1:n_aia)
584 case(335)
585 n_table=n_aia
586 allocate(t_table(1:n_table))
587 allocate(f_table(1:n_table))
588 t_table(1:n_table)=t_aia(1:n_aia)
589 f_table(1:n_table)=f_335(1:n_aia)
590 case(1354)
591 n_table=n_iris
592 allocate(t_table(1:n_table))
593 allocate(f_table(1:n_table))
594 t_table(1:n_table)=t_iris(1:n_iris)
595 f_table(1:n_table)=f_1354(1:n_iris)
596 case(263)
597 n_table=n_eis
598 allocate(t_table(1:n_table))
599 allocate(f_table(1:n_table))
600 t_table(1:n_table)=t_eis1(1:n_eis)
601 f_table(1:n_table)=f_263(1:n_eis)
602 case(264)
603 n_table=n_eis
604 allocate(t_table(1:n_table))
605 allocate(f_table(1:n_table))
606 t_table(1:n_table)=t_eis2(1:n_eis)
607 f_table(1:n_table)=f_264(1:n_eis)
608 case(192)
609 n_table=n_eis
610 allocate(t_table(1:n_table))
611 allocate(f_table(1:n_table))
612 t_table(1:n_table)=t_eis2(1:n_eis)
613 f_table(1:n_table)=f_192(1:n_eis)
614 case(255)
615 n_table=n_eis
616 allocate(t_table(1:n_table))
617 allocate(f_table(1:n_table))
618 t_table(1:n_table)=t_eis2(1:n_eis)
619 f_table(1:n_table)=f_255(1:n_eis)
620 case default
621 allocate(t_table(1))
622 allocate(f_table(1))
623 call mpistop("Unknown wavelength")
624 end select
625 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
626 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
627 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
628 te(ixo^s)=pth(ixo^s)/(ne(ixo^s)*te(ixo^s))*unit_temperature
629 if (si_unit) then
630 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6 ! m^-3 -> cm-3
631 flux(ixo^s)=ne(ixo^s)**2
632 else
633 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
634 flux(ixo^s)=ne(ixo^s)**2
635 endif
636
637 select case(wl)
638 case(94,131,171,193,211,304,335,1354)
639 ! temperature table in log
640 do i=1,n_table
641 if(f_table(i) .gt. 1.d-99) then
642 f_table(i)=log10(f_table(i))
643 else
644 f_table(i)=-99.d0
645 endif
646 enddo
647 loggt=zero
648 {do ix^db=ixomin^db,ixomax^db\}
649 logt=log10(te(ix^d))
650 if (logt>=t_table(1) .and. logt<=t_table(n_table)) then
651 do itt=1,n_table-1
652 if (logt>=t_table(itt) .and. logt<t_table(itt+1)) then
653 loggt=f_table(itt)*(logt-t_table(itt+1))/(t_table(itt)-t_table(itt+1))+&
654 f_table(itt+1)*(logt-t_table(itt))/(t_table(itt+1)-t_table(itt))
655 endif
656 enddo
657 flux(ix^d)=flux(ix^d)*(10**(loggt))
658 if(flux(ix^d)<smalldouble) flux(ix^d)=0.d0
659 else
660 flux(ix^d)=zero
661 endif
662 {enddo\}
663 case default
664 ! temperature table linear
665 {do ix^db=ixomin^db,ixomax^db\}
666 if (te(ix^d)>=t_table(1) .and. te(ix^d)<=t_table(n_table)) then
667 do itt=1,n_table-1
668 if (te(ix^d)>=t_table(itt) .and. te(ix^d)<t_table(itt+1)) then
669 gt=f_table(itt)*(te(ix^d)-t_table(itt+1))/(t_table(itt)-t_table(itt+1))+&
670 f_table(itt+1)*(te(ix^d)-t_table(itt))/(t_table(itt+1)-t_table(itt))
671 endif
672 enddo
673 flux(ix^d)=flux(ix^d)*gt
674 if(flux(ix^d)<smalldouble) flux(ix^d)=0.d0
675 else
676 flux(ix^d)=zero
677 endif
678 {enddo\}
679 end select
680
681 deallocate(t_table,f_table)
682 end subroutine get_euv
683
684 subroutine get_sxr(ixI^L,ixO^L,w,x,fl,flux,El,Eu)
685 !synthesize thermal SXR from El keV to Eu keV released by cm^-3/m^-3
686 ! volume of plasma in 1 s
687 !flux (cgs): photons cm^-3 s^-1
688 !flux (SI): photons m^-3 s^-1
690
691 integer, intent(in) :: ixI^L,ixO^L
692 integer, intent(in) :: El,Eu
693 double precision, intent(in) :: x(ixI^S,1:ndim)
694 double precision, intent(in) :: w(ixI^S,nw)
695 type(te_fluid), intent(in) :: fl
696 double precision, intent(out) :: flux(ixI^S)
697
698 integer :: ix^D,ixO^D
699 integer :: iE,numE
700 double precision :: I0,kb,keV,dE,Ei
701 double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
702 double precision :: Ne(ixI^S),gff(ixI^S),fi(ixI^S)
703 double precision :: EM(ixI^S)
704
705 i0=3.01d-15 ! I0*4*pi*AU**2, I0 from Pinto (2015)
706 kb=const_kb
707 kev=1.0d3*const_ev
708 de=0.1
709 nume=floor((eu-el)/de)
710 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
711 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
712 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
713 te(ixo^s)=pth(ixo^s)/(ne(ixo^s)*te(ixo^s))*unit_temperature
714 if (si_unit) then
715 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6 ! m^-3 -> cm-3
716 em(ixo^s)=(ne(ixo^s))**2*1.d6 ! cm^-3 m^-3
717 else
718 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
719 em(ixo^s)=(ne(ixo^s))**2
720 endif
721 kbt(ixo^s)=kb*te(ixo^s)/kev
722 flux(ixo^s)=0.0d0
723 do ie=0,nume-1
724 ei=de*ie+el*1.d0
725 gff(ixo^s)=1.d0
726 {do ix^db=ixomin^db,ixomax^db\}
727 if (kbt(ix^d)>0.01*ei) then
728 if(kbt(ix^d)<ei) gff(ix^d)=(kbt(ix^d)/ei)**0.4
729 fi(ix^d)=(em(ix^d)*gff(ix^d))*dexp(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
730 else
731 fi(ix^d)=zero
732 endif
733 {enddo\}
734 flux(ixo^s)=flux(ixo^s)+fi(ixo^s)*de
735 enddo
736 flux(ixo^s)=flux(ixo^s)*i0
737 end subroutine get_sxr
738
739 subroutine get_goes_sxr_flux(xbox^L,fl,eflux)
740 !get GOES SXR 1-8A flux observing at 1AU from given box [w/m^2]
742
743 double precision, intent(in) :: xbox^L
744 type(te_fluid), intent(in) :: fl
745 double precision, intent(out) :: eflux
746
747 double precision :: dxb^D,xb^L
748 integer :: iigrid,igrid,j
749 integer :: ixO^L,ixI^L,ix^D
750 double precision :: eflux_grid,eflux_pe
751
752 ^d&iximin^d=ixglo^d;
753 ^d&iximax^d=ixghi^d;
754 ^d&ixomin^d=ixmlo^d;
755 ^d&ixomax^d=ixmhi^d;
756 eflux_pe=zero
757 do iigrid=1,igridstail; igrid=igrids(iigrid);
758 ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
759 ^d&xbmin^d=rnode(rpxmin^d_,igrid);
760 ^d&xbmax^d=rnode(rpxmax^d_,igrid);
761 call get_goes_flux_grid(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,ps(igrid)%dvolume(ixi^s),xbox^l,xb^l,fl,eflux_grid)
762 eflux_pe=eflux_pe+eflux_grid
763 enddo
764 call mpi_allreduce(eflux_pe,eflux,1,mpi_double_precision,mpi_sum,icomm,ierrmpi)
765 end subroutine get_goes_sxr_flux
766
767 subroutine get_goes_flux_grid(ixI^L,ixO^L,w,x,dV,xbox^L,xb^L,fl,eflux_grid)
769
770 integer, intent(in) :: ixI^L,ixO^L
771 double precision, intent(in) :: x(ixI^S,1:ndim),dV(ixI^S)
772 double precision, intent(in) :: w(ixI^S,nw)
773 double precision, intent(in) :: xbox^L,xb^L
774 type(te_fluid), intent(in) :: fl
775 double precision, intent(out) :: eflux_grid
776
777 integer :: ix^D,ixO^D,ixb^L
778 integer :: iE,numE,j,inbox
779 double precision :: I0,kb,keV,dE,Ei,El,Eu,A_cgs
780 double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
781 double precision :: Ne(ixI^S),EM(ixI^S)
782 double precision :: gff,fi,erg_SI
783
784 ! check whether the grid is inside given box
785 inbox=0
786 {if (xbmin^d<xboxmax^d .and. xbmax^d>xboxmin^d) inbox=inbox+1\}
787
788 if (inbox==ndim) then
789 ! indexes for cells inside given box
790 ^d&ixbmin^d=ixomin^d;
791 ^d&ixbmax^d=ixomax^d;
792 {if (xbmax^d>xboxmax^d) ixbmax^d=ixomax^d-ceiling((xbmax^d-xboxmax^d)/dxlevel(^d))\}
793 {if (xbmin^d<xboxmin^d) ixbmin^d=ceiling((xboxmin^d-xbmin^d)/dxlevel(^d))+ixomin^d\}
794
795 i0=1.07d-38 ! photon flux index for observed at 1AU [photon cm^3 m^-2 s^-1 keV^-1]
796 kb=const_kb
797 kev=1.0d3*const_ev
798 erg_si=1.d-7
799 a_cgs=1.d-8 ! Angstrom
800 el=const_h*const_c/(8.d0*a_cgs)/kev ! 8 A
801 eu=const_h*const_c/(1.d0*a_cgs)/kev ! 1 A
802 de=0.1 ! keV
803 nume=floor((eu-el)/de)
804 call fl%get_pthermal(w,x,ixi^l,ixb^l,pth)
805 call fl%get_rho(w,x,ixi^l,ixb^l,ne)
806 call fl%get_var_Rfactor(w,x,ixi^l,ixb^l,te)
807 te(ixb^s)=pth(ixb^s)/(ne(ixb^s)*te(ixb^s))*unit_temperature
808 if (si_unit) then
809 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6 ! m^-3 -> cm-3
810 em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*(unit_length*1.d2)**3 ! cm^-3
811 else
812 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
813 em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*unit_length**3
814 endif
815 kbt(ixb^s)=kb*te(ixb^s)/kev
816 eflux_grid=0.0d0
817
818 do ie=0,nume-1
819 ei=de*ie+el
820 {do ix^db=ixbmin^db,ixbmax^db\}
821 if (kbt(ix^d)>1.d-2*ei) then
822 if(kbt(ix^d)<ei) then
823 gff=(kbt(ix^d)/ei)**0.4
824 else
825 gff=1.d0
826 endif
827 fi=(em(ix^d)*gff)*dexp(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
828 eflux_grid=eflux_grid+fi*de*ei
829 endif
830 {enddo\}
831 enddo
832 eflux_grid=eflux_grid*kev*erg_si
833 endif
834
835 end subroutine get_goes_flux_grid
836
837 {^ifthreed
838 subroutine get_euv_spectrum(qunit,fl)
840
841 integer, intent(in) :: qunit
842 type(te_fluid), intent(in) :: fl
843 character(20) :: datatype
844
845 integer :: mass
846 character (30) :: ion
847 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
848 double precision :: xslit,arcsec
849
850 datatype='spectrum_euv'
851 arcsec=7.25d7/unit_length
852 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
853
854 if (mype==0) print *, '###################################################'
855 select case(spectrum_wl)
856 case (1354)
857 if (mype==0) print *, 'Systhesizing EUV spectrum (observed by IRIS).'
858 case (263,264,192,255)
859 if (mype==0) print *, 'Systhesizing EUV spectrum (observed by Hinode/EIS).'
860 case default
861 call mpistop('Wrong wavelength!')
862 end select
863
865 call mpistop('Wrong spectrum window!')
866 endif
867
868 if (mype==0) write(*,'(a,f8.3,a)') ' Wavelength: ',linecent,' Angstrom'
869 if (mype==0) print *, 'Unit of EUV flux: DN s^-1 pixel^-1'
870
871 if (dat_resolution) then
872 if (mype==0) then
873 write(*,'(a,f5.3,a,f5.1,a)') ' Supposed pixel: ',wlrsl,' Angstrom x ',spacersl*725.0, ' km'
874 print *, 'Unit of wavelength: Angstrom (0.1 nm) '
875 if (si_unit) then
876 write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d6,' Mm'
877 else
878 write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d8,' Mm'
879 endif
880 write(*,'(a,f8.1,a)') ' Supposed width of slit: ',wslit*725.0,' km'
881 endif
882 call get_spectrum_datresol(qunit,datatype,fl)
883 else
884 if (mype==0) then
885 print *, 'Unit of wavelength: Angstrom (0.1 nm) '
886 if (activate_unit_arcsec) then
887 write(*,'(a,f5.3,a,f5.1,a)') ' Pixel: ',wlrsl,' Angstrom x ',spacersl*725.0, ' km'
888 print *, 'Unit of length: arcsec (~725 km)'
889 write(*,'(a,f8.1,a)') ' Location of slit: xI1 = ',location_slit,' arcsec'
890 write(*,'(a,f8.1,a)') ' Width of slit: ',wslit,' arcsec'
891 else
892 if (si_unit) then
893 if (mype==0) write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d6,' Mm'
894 else
895 if (mype==0) write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d8,' Mm'
896 endif
897 write(*,'(a,f8.1,a)') ' Location of slit: xI1 = ',location_slit,' Unit_length'
898 write(*,'(a,f8.1,a)') ' Width of slit: ',wslit*725.d0,' km'
899 endif
900 endif
901 if (mype==0) print *, 'Direction of the slit: parallel to xI2 vector'
903 call get_spectrum(qunit,datatype,fl)
904 else
905 call mpistop("EUV spectrum synthesis: support for sperical coordinates is to be added!")
906 endif
907 endif
908
909 if (mype==0) print *, '###################################################'
910
911 end subroutine get_euv_spectrum
912
913 subroutine get_spectrum_datresol(qunit,datatype,fl)
914
915 integer, intent(in) :: qunit
916 character(20), intent(in) :: datatype
917 type(te_fluid), intent(in) :: fl
918
919 integer :: numWL,numXS,iwL,ixS,numWI,numS
920 double precision :: dwLg,xSmin,xSmax,wLmin,wLmax
921 double precision, allocatable :: wL(:),xS(:),dwL(:),dxS(:)
922 double precision, allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
923 integer :: strtype,nstrb,nbb,nuni,nstr,bnx
924 double precision :: qs,dxfirst,dxmid,lenstr
925
926 integer :: iigrid,igrid,j,dir_loc
927 double precision :: xbmin(1:ndim),xbmax(1:ndim)
928
929 dwlg=1.d-3
930 numwl=4*int((spectrum_window_max-spectrum_window_min)/(4.d0*dwlg))
931 wlmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwlg*numwl/2
932 wlmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwlg*numwl/2
933 allocate(wl(numwl),dwl(numwl))
934 dwl(:)=dwlg
935 do iwl=1,numwl
936 wl(iwl)=wlmin+iwl*dwlg-half*dwlg
937 enddo
938
939 select case(direction_slit)
940 case (1)
941 numxs=domain_nx1*2**(refine_max_level-1)
942 xsmin=xprobmin1
943 xsmax=xprobmax1
944 bnx=block_nx1
945 nbb=domain_nx1
946 strtype=stretch_type(1)
947 nstrb=nstretchedblocks_baselevel(1)
948 qs=qstretch_baselevel(1)
949 if (mype==0) print *, 'Direction of the slit: x'
950 case (2)
951 numxs=domain_nx2*2**(refine_max_level-1)
952 xsmin=xprobmin2
953 xsmax=xprobmax2
954 bnx=block_nx2
955 nbb=domain_nx2
956 strtype=stretch_type(2)
957 nstrb=nstretchedblocks_baselevel(2)
958 qs=qstretch_baselevel(2)
959 if (mype==0) print *, 'Direction of the slit: y'
960 case (3)
961 numxs=domain_nx3*2**(refine_max_level-1)
962 xsmin=xprobmin3
963 xsmax=xprobmax3
964 bnx=block_nx3
965 nbb=domain_nx3
966 strtype=stretch_type(3)
967 nstrb=nstretchedblocks_baselevel(3)
968 qs=qstretch_baselevel(3)
969 if (mype==0) print *, 'Direction of the slit: z'
970 case default
971 call mpistop('Wrong direction_slit')
972 end select
973
974 allocate(xs(numxs),dxs(numxs),spectra(numwl,numxs),spectra_rc(numwl,numxs))
975 numwi=1
976 allocate(wi(numwl,numxs,numwi))
977
978 select case(strtype)
979 case(0) ! uniform
980 dxs(:)=(xsmax-xsmin)/numxs
981 do ixs=1,numxs
982 xs(ixs)=xsmin+dxs(ixs)*(ixs-half)
983 enddo
984 case(1) ! uni stretch
985 qs=qs**(one/2**(refine_max_level-1))
986 dxfirst=(xsmax-xsmin)*(one-qs)/(one-qs**numxs)
987 dxs(1)=dxfirst
988 do ixs=2,numxs
989 dxs(ixs)=dxfirst*qs**(ixs-1)
990 xs(ixs)=dxs(1)/(one-qs)*(one-qs**(ixs-1))+half*dxs(ixs)
991 enddo
992 case(2) ! symm stretch
993 ! base level, nbb = nstr + nuni + nstr
994 nstr=nstrb*bnx/2
995 nuni=nbb-nstrb*bnx
996 lenstr=(xsmax-xsmin)/(2.d0+nuni*(one-qs)/(one-qs**nstr))
997 dxfirst=(xsmax-xsmin)/(dble(nuni)+2.d0/(one-qs)*(one-qs**nstr))
998 dxmid=dxfirst
999 ! refine_max level, numXI = nstr + nuni + nstr
1000 nstr=nstr*2**(refine_max_level-1)
1001 nuni=nuni*2**(refine_max_level-1)
1002 qs=qs**(one/2**(refine_max_level-1))
1003 dxfirst=lenstr*(one-qs)/(one-qs**nstr)
1004 dxmid=dxmid/2**(refine_max_level-1)
1005 ! uniform center
1006 if(nuni .gt. 0) then
1007 do ixs=nstr+1,nstr+nuni
1008 dxs(ixs)=dxmid
1009 xs(ixs)=lenstr+(dble(ixs)-0.5d0-nstr)*dxs(ixs)+xsmin
1010 enddo
1011 endif
1012 ! left half
1013 do ixs=nstr,1,-1
1014 dxs(ixs)=dxfirst*qs**(nstr-ixs)
1015 xs(ixs)=xsmin+lenstr-dxs(ixs)*half-dxfirst*(one-qs**(nstr-ixs))/(one-qs)
1016 enddo
1017 ! right half
1018 do ixs=nstr+nuni+1,numxs
1019 dxs(ixs)=dxfirst*qs**(ixs-nstr-nuni-1)
1020 xs(ixs)=xsmax-lenstr+dxs(ixs)*half+dxfirst*(one-qs**(ixs-nstr-nuni-1))/(one-qs)
1021 enddo
1022 case default
1023 call mpistop("unknown stretch type")
1024 end select
1025
1026 if (los_phi==0 .and. los_theta==90 .and. direction_slit==2) then
1027 ! LOS->x slit->y
1028 dir_loc=3
1029 else if (los_phi==0 .and. los_theta==90 .and. direction_slit==3) then
1030 ! LOS->x slit->z
1031 dir_loc=2
1032 else if (los_phi==90 .and. los_theta==90 .and. direction_slit==1) then
1033 ! LOS->y slit->x
1034 dir_loc=3
1035 else if (los_phi==90 .and. los_theta==90 .and. direction_slit==3) then
1036 ! LOS->y slit->z
1037 dir_loc=1
1038 else if (los_theta==0 .and. direction_slit==1) then
1039 ! LOS->z slit->x
1040 dir_loc=2
1041 else if (los_theta==0 .and. direction_slit==2) then
1042 ! LOS->z slit->y
1043 dir_loc=1
1044 else
1045 call mpistop('Wrong combination of LOS and slit direction!')
1046 endif
1047
1048 if (dir_loc==1) then
1049 if (location_slit>xprobmax1 .or. location_slit<xprobmin1) then
1050 call mpistop('Wrong value for location_slit!')
1051 endif
1052 if(mype==0) write(*,'(a,f8.1,a)') ' Location of slit: x = ',location_slit,' Unit_length'
1053 else if (dir_loc==2) then
1054 if (location_slit>xprobmax2 .or. location_slit<xprobmin2) then
1055 call mpistop('Wrong value for location_slit!')
1056 endif
1057 if(mype==0) write(*,'(a,f8.1,a)') ' Location of slit: y = ',location_slit,' Unit_length'
1058 else
1059 if (location_slit>xprobmax3 .or. location_slit<xprobmin3) then
1060 call mpistop('Wrong value for location_slit!')
1061 endif
1062 if(mype==0) write(*,'(a,f8.1,a)') ' Location of slit: z = ',location_slit,' Unit_length'
1063 endif
1064
1065 ! find slit and do integration
1066 spectra=zero
1067 do iigrid=1,igridstail; igrid=igrids(iigrid);
1068 ^d&xbmin(^d)=rnode(rpxmin^d_,igrid);
1069 ^d&xbmax(^d)=rnode(rpxmax^d_,igrid);
1070 if (location_slit>=xbmin(dir_loc) .and. location_slit<xbmax(dir_loc)) then
1071 call integrate_spectra_datresol(igrid,wl,dwl,spectra,numwl,numxs,dir_loc,fl)
1072 endif
1073 enddo
1074
1075 nums=numwl*numxs
1076 call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
1077 mpi_sum,icomm,ierrmpi)
1078 do iwl=1,numwl
1079 do ixs=1,numxs
1080 if (spectra_rc(iwl,ixs)>smalldouble) then
1081 wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
1082 else
1083 wi(iwl,ixs,1)=zero
1084 endif
1085 enddo
1086 enddo
1087
1088 call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
1089
1090 deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
1091
1092 end subroutine get_spectrum_datresol
1093
1094 subroutine integrate_spectra_datresol(igrid,wL,dwL,spectra,numWL,numXS,dir_loc,fl)
1095 use mod_constants
1096
1097 integer, intent(in) :: igrid,numWL,numXS,dir_loc
1098 type(te_fluid), intent(in) :: fl
1099 double precision, intent(in) :: wL(numWL),dwL(numWL)
1100 double precision, intent(inout) :: spectra(numWL,numXS)
1101
1102 integer :: direction_LOS
1103 integer :: ixO^L,ixI^L,ix^D,ixOnew
1104 double precision, allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
1105 double precision :: wlc,wlwd
1106
1107 integer :: mass
1108 double precision :: logTe,lineCent
1109 character (30) :: ion
1110 double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
1111
1112 integer :: levelg,rft,ixSmin,ixSmax,iwL
1113 double precision :: flux_pix,dL
1114
1115 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1116
1117 if (los_phi==0 .and. los_theta==90) then
1118 direction_los=1
1119 else if (los_phi==90 .and. los_theta==90) then
1120 direction_los=2
1121 else
1122 direction_los=3
1123 endif
1124
1125 ^d&ixomin^d=ixmlo^d\
1126 ^d&ixomax^d=ixmhi^d\
1127 ^d&iximin^d=ixglo^d\
1128 ^d&iximax^d=ixghi^d\
1129 allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
1130
1131 ^d&ix^d=ixomin^d;
1132 if (dir_loc==1) then
1133 do ix1=ixomin1,ixomax1
1134 if (location_slit>=(ps(igrid)%x(ix^d,1)-half*ps(igrid)%dx(ix^d,1)) .and. &
1135 location_slit<(ps(igrid)%x(ix^d,1)+half*ps(igrid)%dx(ix^d,1))) then
1136 ixonew=ix1
1137 endif
1138 enddo
1139 ixomin1=ixonew
1140 ixomax1=ixonew
1141 else if (dir_loc==2) then
1142 do ix2=ixomin2,ixomax2
1143 if (location_slit>=(ps(igrid)%x(ix^d,2)-half*ps(igrid)%dx(ix^d,2)) .and. &
1144 location_slit<(ps(igrid)%x(ix^d,2)+half*ps(igrid)%dx(ix^d,2))) then
1145 ixonew=ix2
1146 endif
1147 enddo
1148 ixomin2=ixonew
1149 ixomax2=ixonew
1150 else
1151 do ix3=ixomin3,ixomax3
1152 if (location_slit>=(ps(igrid)%x(ix^d,3)-half*ps(igrid)%dx(ix^d,3)) .and. &
1153 location_slit<(ps(igrid)%x(ix^d,3)+half*ps(igrid)%dx(ix^d,3))) then
1154 ixonew=ix3
1155 endif
1156 enddo
1157 ixomin3=ixonew
1158 ixomax3=ixonew
1159 endif
1160
1161 call get_euv(spectrum_wl,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
1162 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2 ! adjust flux due to artifical change of resolution
1163 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
1164 v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
1165 call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
1166 call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
1167 te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
1168
1169 ! grid parameters
1170 levelg=ps(igrid)%level
1171 rft=2**(refine_max_level-levelg)
1172
1173 {do ix^d=ixomin^d,ixomax^d\}
1174 if (flux(ix^d)>smalldouble) then
1175 if (si_unit) then
1176 wlc=linecent*(1.d0+v(ix^d)*unit_velocity*1.d2/const_c)
1177 else
1178 wlc=linecent*(1.d0+v(ix^d)*unit_velocity/const_c)
1179 endif
1180 wlwd=sqrt(kb_cgs*te(ix^d)*unit_temperature/(mass*mp_cgs))
1181 wlwd=wlwd*linecent/const_c
1182 ! involved pixel
1183 select case(direction_slit)
1184 case(1)
1185 ixsmin=(block_nx1*(node(pig1_,igrid)-1)+(ix1-ixomin1))*rft+1
1186 ixsmax=(block_nx1*(node(pig1_,igrid)-1)+(ix1-ixomin1+1))*rft
1187 case(2)
1188 ixsmin=(block_nx2*(node(pig2_,igrid)-1)+(ix2-ixomin2))*rft+1
1189 ixsmax=(block_nx2*(node(pig2_,igrid)-1)+(ix2-ixomin2+1))*rft
1190 case(3)
1191 ixsmin=(block_nx3*(node(pig3_,igrid)-1)+(ix3-ixomin3))*rft+1
1192 ixsmax=(block_nx3*(node(pig3_,igrid)-1)+(ix3-ixomin3+1))*rft
1193 end select
1194 ! LOS depth
1195 select case(direction_los)
1196 case(1)
1197 dl=ps(igrid)%dx(ix^d,1)*unit_length
1198 case(2)
1199 dl=ps(igrid)%dx(ix^d,2)*unit_length
1200 case default
1201 dl=ps(igrid)%dx(ix^d,3)*unit_length
1202 end select
1203 if (si_unit) dl=dl*1.d2
1204 ! integral pixel flux
1205 do iwl=1,numwl
1206 flux_pix=flux(ix^d)*wlrsl*dl*exp(-(wl(iwl)-wlc)**2/(2*wlwd**2))/(sqrt(2*dpi)*wlwd)
1207 if (flux_pix>smalldouble) then
1208 flux_pix=flux_pix*wslit/spacersl
1209 spectra(iwl,ixsmin:ixsmax)=spectra(iwl,ixsmin:ixsmax)+flux_pix
1210 endif
1211 enddo
1212 endif
1213 {enddo\}
1214
1215 deallocate(flux,v,pth,te,rho)
1216
1217 end subroutine integrate_spectra_datresol
1218
1219 subroutine get_spectrum(qunit,datatype,fl)
1220
1221 integer, intent(in) :: qunit
1222 character(20), intent(in) :: datatype
1223 type(te_fluid), intent(in) :: fl
1224
1225 integer :: numWL,numXS,iwL,ixS,numWI,ix^D
1226 double precision :: dwLg,dxSg,xSmin,xSmax,xScent,wLmin,wLmax
1227 double precision, allocatable :: wL(:),xS(:),dwL(:),dxS(:)
1228 double precision, allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
1229 double precision :: vec_cor(1:3),xI_cor(1:2)
1230 double precision :: res,r_loc,r_max
1231
1232 integer :: mass
1233 character (30) :: ion
1234 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1235 double precision :: unitv,arcsec,RHESSI_rsl,pixel
1236 integer :: iigrid,igrid,i,j,numS
1237 double precision :: xLmin,xLmax,xslit
1238
1239 if (coordinate==spherical) then
1241 else
1242 ! cartesian
1244 endif
1245
1246 ! calculate domain in space
1247 if (coordinate==spherical) then
1248 xsmin=-abs(xprobmax1)
1249 xsmax=abs(xprobmax1)
1250 else
1251 do ix1=1,2
1252 if (ix1==1) vec_cor(1)=xprobmin1
1253 if (ix1==2) vec_cor(1)=xprobmax1
1254 do ix2=1,2
1255 if (ix2==1) vec_cor(2)=xprobmin2
1256 if (ix2==2) vec_cor(2)=xprobmax2
1257 do ix3=1,2
1258 if (ix3==1) vec_cor(3)=xprobmin3
1259 if (ix3==2) vec_cor(3)=xprobmax3
1260 if (big_image) then
1261 r_loc=(vec_cor(1)-x_origin(1))**2
1262 r_loc=r_loc+(vec_cor(2)-x_origin(2))**2
1263 r_loc=r_loc+(vec_cor(3)-x_origin(3))**2
1264 r_loc=sqrt(r_loc)
1265 if (ix1==1 .and. ix2==1 .and. ix3==1) then
1266 r_max=r_loc
1267 else
1268 r_max=max(r_max,r_loc)
1269 endif
1270 else
1271 call get_cor_image(vec_cor,xi_cor)
1272 if (ix1==1 .and. ix2==1 .and. ix3==1) then
1273 xsmin=xi_cor(2)
1274 xsmax=xi_cor(2)
1275 else
1276 xsmin=min(xsmin,xi_cor(2))
1277 xsmax=max(xsmax,xi_cor(2))
1278 endif
1279 endif
1280 enddo
1281 enddo
1282 enddo
1283 if (big_image) then
1284 xsmin=-r_max
1285 xsmax=r_max
1286 endif
1287 endif
1288 xscent=(xsmin+xsmax)/2.d0
1289
1290 ! tables for storing spectra data
1291 if (si_unit) then
1292 arcsec=7.25d5/unit_length
1293 else
1294 arcsec=7.25d7/unit_length
1295 endif
1296 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1297 dxsg=spacersl*arcsec
1298 numxs=ceiling((xsmax-xscent)/dxsg)
1299 xsmin=xscent-numxs*dxsg
1300 xsmax=xscent+numxs*dxsg
1301 numxs=numxs*2
1302 dwlg=wlrsl
1303 numwl=2*int((spectrum_window_max-spectrum_window_min)/(2.d0*dwlg))
1304 wlmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwlg*numwl/2
1305 wlmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwlg*numwl/2
1306 allocate(wl(numwl),dwl(numwl),xs(numxs),dxs(numxs))
1307 numwi=1
1308 allocate(wi(numwl,numxs,numwi),spectra(numwl,numxs),spectra_rc(numwl,numxs))
1309 do iwl=1,numwl
1310 wl(iwl)=wlmin+iwl*dwlg-half*dwlg
1311 dwl=dwlg
1312 enddo
1313 do ixs=1,numxs
1314 xs(ixs)=xsmin+dxsg*(ixs-half)
1315 dxs(ixs)=dxsg
1316 enddo
1317
1318 ! find slit and do integration
1319 spectra=zero
1320 do iigrid=1,igridstail; igrid=igrids(iigrid);
1321 do ix1=1,2
1322 if (ix1==1) vec_cor(1)=rnode(rpxmin1_,igrid)
1323 if (ix1==2) vec_cor(1)=rnode(rpxmax1_,igrid)
1324 do ix2=1,2
1325 if (ix2==1) vec_cor(2)=rnode(rpxmin2_,igrid)
1326 if (ix2==2) vec_cor(2)=rnode(rpxmax2_,igrid)
1327 do ix3=1,2
1328 if (ix3==1) vec_cor(3)=rnode(rpxmin3_,igrid)
1329 if (ix3==2) vec_cor(3)=rnode(rpxmax3_,igrid)
1330 call get_cor_image(vec_cor,xi_cor)
1331 if (ix1==1 .and. ix2==1 .and. ix3==1) then
1332 xlmin=xi_cor(1)
1333 xlmax=xi_cor(1)
1334 else
1335 xlmin=min(xlmin,xi_cor(1))
1336 xlmax=max(xlmax,xi_cor(1))
1337 endif
1338 enddo
1339 enddo
1340 enddo
1341
1342 if (activate_unit_arcsec) then
1343 xslit=location_slit*arcsec
1344 else
1345 xslit=location_slit
1346 endif
1347 if (xslit>=xlmin-wslit*arcsec .and. xslit<=xlmax+wslit*arcsec) then
1348 call integrate_spectra_cartesian(igrid,wl,dwlg,xs,dxsg,spectra,numwl,numxs,fl)
1349 endif
1350 enddo
1351
1352 nums=numwl*numxs
1353 call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
1354 mpi_sum,icomm,ierrmpi)
1355 do iwl=1,numwl
1356 do ixs=1,numxs
1357 if (spectra_rc(iwl,ixs)>smalldouble) then
1358 wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
1359 else
1360 wi(iwl,ixs,1)=zero
1361 endif
1362 enddo
1363 enddo
1364
1365 if (activate_unit_arcsec) then
1366 xs=xs/arcsec
1367 dxs=dxs/arcsec
1368 endif
1369
1370 call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
1371
1372 deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
1373
1374 end subroutine get_spectrum
1375
1376 subroutine integrate_spectra_cartesian(igrid,wL,dwLg,xS,dxSg,spectra,numWL,numXS,fl)
1377
1378 integer, intent(in) :: igrid,numWL,numXS
1379 double precision, intent(in) :: wL(numWL),xS(numXS)
1380 double precision, intent(in) :: dwLg,dxSg
1381 double precision, intent(inout) :: spectra(numWL,numXS)
1382 type(te_fluid), intent(in) :: fl
1383
1384 integer :: ixO^L,ixI^L,ix^D,ixOnew,j
1385 double precision, allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
1386 double precision :: wlc,wlwd,res,dst_slit,xslit,arcsec
1387 double precision :: vloc(1:3),xloc(1:3),dxloc(1:3),xIloc(1:2),dxIloc(1:2)
1388 integer :: nSubC^D,iSubC^D,iwL,ixS,ixSmin,ixSmax,iwLmin,iwLmax,nwL
1389 double precision :: slit_width,dxSubC^D,xerf^L,fluxSubC
1390 double precision :: xSubC(1:3),xCent(1:2)
1391
1392 integer :: mass
1393 double precision :: logTe,lineCent
1394 character (30) :: ion
1395 double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
1396 double precision :: sigma_wl,sigma_xs,factor
1397
1398 if (si_unit) then
1399 arcsec=7.25d5/unit_length
1400 else
1401 arcsec=7.25d7/unit_length
1402 endif
1403 if (activate_unit_arcsec) then
1404 xslit=location_slit*arcsec
1405 else
1406 xslit=location_slit
1407 endif
1408
1409 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1410
1411 ^d&ixomin^d=ixmlo^d\
1412 ^d&ixomax^d=ixmhi^d\
1413 ^d&iximin^d=ixglo^d\
1414 ^d&iximax^d=ixghi^d\
1415 allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
1416 ! get local EUV flux and velocity
1417 call get_euv(spectrum_wl,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
1418 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2 ! adjust flux due to artifical change of resolution
1419 call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
1420 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
1421 call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
1422 te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
1423 {do ix^d=ixomin^d,ixomax^d\}
1424 do j=1,3
1425 vloc(j)=ps(igrid)%w(ix^d,iw_mom(j))/rho(ix^d)
1426 enddo
1427 call dot_product_loc(vloc,vec_los,res)
1428 v(ix^d)=res
1429 {enddo\}
1430
1431 deallocate(rho)
1432
1433 slit_width=wslit*arcsec
1434 sigma_wl=sigma_psf*dwlg
1435 sigma_xs=sigma_psf*dxsg
1436 {do ix^d=ixomin^d,ixomax^d\}
1437 if (flux(ix^d)>smalldouble) then
1438 xloc(1:3)=ps(igrid)%x(ix^d,1:3)
1439 dxloc(1:3)=ps(igrid)%dx(ix^d,1:3)
1440 call get_cor_image(xloc,xiloc)
1441 call dot_product_loc(dxloc,vec_xi1,res)
1442 dxiloc(1)=abs(res)
1443 if (xiloc(1)>=xslit-half*(slit_width+dxiloc(1)) .and. &
1444 xiloc(1)<=xslit+half*(slit_width+dxiloc(1))) then
1445 ^d&nsubc^d=1;
1446 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(vec_xi1(^d))/(slit_width/16.d0)));
1447 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(vec_xi2(^d))/(dxsg/4.d0)));
1448 ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
1449 ! local line center and line width
1450 if (si_unit) then
1451 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length*1.d2/dxsg/dxsg ! DN s^-1
1452 wlc=linecent*(1.d0+v(ix^d)*unit_velocity*1.d2/const_c)
1453 else
1454 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length/dxsg/dxsg ! DN s^-1
1455 wlc=linecent*(1.d0+v(ix^d)*unit_velocity/const_c)
1456 endif
1457 wlwd=sqrt(kb_cgs*te(ix^d)*unit_temperature/(mass*mp_cgs))
1458 wlwd=wlwd*linecent/const_c
1459 ! dividing a cell to several parts to get more accurate integrating values
1460 {do isubc^d=1,nsubc^d\}
1461 ^d&xsubc(^d)=xloc(^d)-half*dxloc(^d)+(isubc^d-half)*dxsubc^d;
1462 call get_cor_image(xsubc,xcent)
1463 dst_slit=abs(xcent(1)-xslit) ! space distance to slit center
1464 if (dst_slit<=half*slit_width) then
1465 ixs=floor((xcent(2)-(xs(1)-half*dxsg))/dxsg)+1
1466 ixsmin=max(1,ixs-3)
1467 ixsmax=min(ixs+3,numxs)
1468 iwl=floor((wlc-(wl(1)-half*dwlg))/dwlg)+1
1469 nwl=3*ceiling(wlwd/dwlg+1)
1470 iwlmin=max(1,iwl-nwl)
1471 iwlmax=min(iwl+nwl,numwl)
1472 ! calculate the contribution to nearby pixels
1473 do iwl=iwlmin,iwlmax
1474 do ixs=ixsmin,ixsmax
1475 xerfmin1=(wl(iwl)-half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
1476 xerfmax1=(wl(iwl)+half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
1477 xerfmin2=(xs(ixs)-half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
1478 xerfmax2=(xs(ixs)+half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
1479 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
1480 spectra(iwl,ixs)=spectra(iwl,ixs)+fluxsubc*factor
1481 enddo
1482 enddo
1483 ! nearby pixels
1484 endif
1485 {enddo\}
1486 endif
1487 endif
1488 {enddo\}
1489
1490 deallocate(flux,v,pth,te)
1491 end subroutine integrate_spectra_cartesian
1492 }
1493
1494 {^ifthreed
1495 subroutine get_euv_image(qunit,fl)
1497
1498 integer, intent(in) :: qunit
1499 type(te_fluid), intent(in) :: fl
1500 character(20) :: datatype
1501
1502 integer :: mass
1503 character (30) :: ion
1504 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1505 double precision :: t0,t1
1506
1507 t0=mpi_wtime()
1508 datatype='image_euv'
1509 call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1510
1511 if (mype==0) then
1512 print *, '###################################################'
1513 print *, 'Systhesizing EUV image'
1514 write(*,'(a,f8.3,a)') ' Wavelength: ',linecent,' Angstrom'
1515 print *, 'Unit of EUV flux: DN s^-1 pixel^-1'
1516 endif
1517
1518 if (dat_resolution) then
1519 if (coordinate/=cartesian) call mpistop('EUV synthesis: only cartesian is supported for .dat resolution!')
1520 if (mype==0) then
1521 write(*,'(a,f7.1,a,f7.1,a,f5.1,a,f5.1,a)') ' Supposed Pixel: ',spacersl*725.0,' km x ',spacersl*725.0, &
1522 ' km (', spacersl, ' arcsec x ', spacersl, ' arcsec)'
1523 if (si_unit) then
1524 write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d6,' Mm'
1525 else
1526 write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d8,' Mm'
1527 endif
1528 endif
1529 if (los_phi==0 .and. los_theta==90) then
1530 call get_image_datresol(qunit,datatype,fl)
1531 else if (los_phi==90 .and. los_theta==90) then
1532 call get_image_datresol(qunit,datatype,fl)
1533 else if (los_theta==0) then
1534 call get_image_datresol(qunit,datatype,fl)
1535 else
1536 call mpistop('ERROR: Wrong LOS for synthesizing emission!')
1537 endif
1538 else
1539 if (mype==0) then
1540 write(*,'(a,f7.1,a,f7.1,a,f5.1,a,f5.1,a)') ' Pixel: ',spacersl*725.0,' km x ',spacersl*725.0, ' km (', &
1541 spacersl, ' arcsec x ', spacersl, ' arcsec)'
1542 if (activate_unit_arcsec) then
1543 print *, 'Unit of length: arcsec (~725 km)'
1544 else
1545 if (si_unit) then
1546 write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d6,' Mm'
1547 else
1548 write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d8,' Mm'
1549 endif
1550 endif
1551 endif
1552 if (coordinate==cartesian) then
1553 if (mype==0) write(*,'(a,f8.3,f8.3,f8.3,a)') ' Mapping: [',x_origin(1),x_origin(2),x_origin(3), &
1554 '] of the simulation box is located at [X=0,Y=0] of the image'
1555 call get_image(qunit,datatype,fl)
1556 else if (coordinate==spherical) then
1557 if (mype==0) write(*,'(a,f6.3,f8.3,f8.3,a)') ' Mapping: R=0 of the simulation box is located at [X=0,Y=0] of the image'
1558 call get_image(qunit,datatype,fl)
1559 else
1560 call mpistop("EUV synthesis: this coordinate is not supported!")
1561 endif
1562 endif
1563
1564 t1=mpi_wtime()
1565 if (mype==0) print *, 'time comsuming: ',t1-t0,' s'
1566 if (mype==0) print *, '###################################################'
1567
1568 end subroutine get_euv_image
1569
1570 subroutine get_sxr_image(qunit,fl)
1572
1573 integer, intent(in) :: qunit
1574 type(te_fluid), intent(in) :: fl
1575 character(20) :: datatype
1576 double precision :: RHESSI_rsl
1577 double precision :: t0,t1
1578
1579 t0=mpi_wtime()
1580 datatype='image_sxr'
1581 rhessi_rsl=2.3/instrument_resolution_factor
1582
1583 if (mype==0) then
1584 print *, '###################################################'
1585 print *, 'Systhesizing SXR image (observed at 1 AU).'
1586 write(*,'(a,i2,a,i2,a)') ' Passband: ',emin_sxr,' - ',emax_sxr,' keV'
1587 endif
1588
1589 if (dat_resolution) then
1590 if (coordinate/=cartesian) call mpistop('SXR synthesis: only cartesian is supported for .dat resolution!')
1591 if (mype==0) then
1592 print *, 'Unit of SXR flux: photons cm^-2 s^-1 pixel^-1'
1593 write(*,'(a,f5.1,a,f5.1,a,f5.1,a,f5.1,a)') ' Supposed Pixel: ',rhessi_rsl*0.725, ' Mm x ',rhessi_rsl*0.725, &
1594 ' Mm (', rhessi_rsl, ' arcsec x ', rhessi_rsl, ' arcsec)'
1595 if (si_unit) then
1596 write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d6,' Mm'
1597 else
1598 write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d8,' Mm'
1599 endif
1600 endif
1601 if (los_phi==0 .and. los_theta==90) then
1602 call get_image_datresol(qunit,datatype,fl)
1603 else if (los_phi==90 .and. los_theta==90) then
1604 call get_image_datresol(qunit,datatype,fl)
1605 else if (los_theta==0) then
1606 call get_image_datresol(qunit,datatype,fl)
1607 else
1608 call mpistop('ERROR: Wrong LOS for synthesizing emission!')
1609 endif
1610 else
1611 if (mype==0) then
1612 print *, 'Unit of SXR flux: photons cm^-2 s^-1 pixel^-1'
1613 write(*,'(a,f5.1,a,f5.1,a,f5.1,a,f5.1,a)') ' Pixel: ',rhessi_rsl*0.725, ' Mm x ',rhessi_rsl*0.725, &
1614 ' Mm (', rhessi_rsl, ' arcsec x ', rhessi_rsl, ' arcsec)'
1615 if (activate_unit_arcsec) then
1616 print *, 'Unit of length: arcsec (~725 km)'
1617 else
1618 if (si_unit) then
1619 write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d6,' Mm'
1620 else
1621 write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d8,' Mm'
1622 endif
1623 endif
1624 endif
1625 if (coordinate==cartesian) then
1626 if (mype==0) write(*,'(a,f8.3,f8.3,f8.3,a)') ' Mapping: [',x_origin(1),x_origin(2),x_origin(3), &
1627 '] of the simulation box is located at [X=0,Y=0] of the image'
1628 call get_image(qunit,datatype,fl)
1629 else if (coordinate==spherical) then
1630 if (mype==0) write(*,'(a,f6.3,f8.3,f8.3,a)') ' Mapping: R=0 of the simulation box is located at [X=0,Y=0] of the image'
1631 call get_image(qunit,datatype,fl)
1632 else
1633 call mpistop("SXR synthesis: this coordinate is not supported!")
1634 endif
1635 endif
1636
1637 t1=mpi_wtime()
1638 if (mype==0) print *, 'time comsuming:',t1-t0
1639 if (mype==0) print *, '###################################################'
1640
1641 end subroutine get_sxr_image
1642
1643 subroutine get_whitelight_image(qunit,fl)
1645
1646 integer, intent(in) :: qunit
1647 type(te_fluid), intent(in) :: fl
1648 character(20) :: datatype
1649 double precision :: LASCO_rsl
1650
1651 if (mype==0) print *, '###################################################'
1652
1653 if (whitelight_instrument=='LASCO/C1') then
1654 lasco_rsl=5.6d0/instrument_resolution_factor
1655 if (mype==0) print *, 'Systhesizing white light image (observed by LASCO/C1).'
1656 else if (whitelight_instrument=='LASCO/C2') then
1657 lasco_rsl=11.4d0/instrument_resolution_factor
1658 if (mype==0) print *, 'Systhesizing white light image (observed by LASCO/C2).'
1659 else if (whitelight_instrument=='LASCO/C3') then
1660 lasco_rsl=56.d0/instrument_resolution_factor
1661 if (mype==0) print *, 'Systhesizing white light image (observed by LASCO/C3).'
1662 else
1663 call mpistop('Whitelight synthesis: instrument is not supported!')
1664 endif
1665
1666 if (mype==0) write(*,'(a,f5.1,a,f5.1,a,f5.1,a,f5.1,a)') ' Pixel: ',lasco_rsl*0.725,' Mm x ',lasco_rsl*0.725, ' Mm (', &
1667 lasco_rsl, ' arcsec x ', lasco_rsl, ' arcsec) '
1668 if (mype==0) print *, 'Unit of white light flux: average Sun brightness'
1669
1670 datatype='image_whitelight'
1671
1672 if (mype==0) then
1673 if (activate_unit_arcsec) then
1674 print *, 'Unit of length: arcsec (~725 km)'
1675 else
1676 if (si_unit) then
1677 if (mype==0) write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d6,' Mm'
1678 else
1679 if (mype==0) write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d8,' Mm'
1680 endif
1681 endif
1682 endif
1683
1684 if (coordinate==spherical) then
1685 if (mype==0) write(*,'(a,f6.3,f8.3,f8.3,a)') ' Mapping: R=0 of the simulation box is located at [X=0,Y=0] of the image'
1686 call get_image(qunit,datatype,fl)
1687 else
1688 call mpistop("Whitelight synthesis: this coordinate is not supported!")
1689 endif
1690
1691 if (mype==0) print *, '###################################################'
1692
1693 end subroutine get_whitelight_image
1694
1695 subroutine get_image_datresol(qunit,datatype,fl)
1696 ! integrate emission flux along line of sight (LOS)
1697 ! in a 3D simulation box and get a 2D EUV image
1699 use mod_constants
1700
1701 integer, intent(in) :: qunit
1702 character(20), intent(in) :: datatype
1703 type(te_fluid), intent(in) :: fl
1704
1705 double precision :: dx^D
1706 integer :: numX^D,ix^D
1707 double precision, allocatable :: EUV(:,:),EUVs(:,:),Dpl(:,:),Dpls(:,:)
1708 double precision, allocatable :: SXR(:,:),SXRs(:,:),wI(:,:,:)
1709 double precision, allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:),dxIi
1710 integer :: numXI1,numXI2,numSI,numWI
1711 double precision :: xI^L
1712 integer :: iigrid,igrid,i,j
1713 double precision, allocatable :: xIF1(:),xIF2(:),dxIF1(:),dxIF2(:)
1714 integer :: nXIF1,nXIF2
1715 double precision :: xIF^L
1716
1717 double precision :: unitv,arcsec,RHESSI_rsl
1718 integer :: strtype^D,nstrb^D,nbb^D,nuni^D,nstr^D,bnx^D
1719 double precision :: qs^D,dxfirst^D,dxmid^D,lenstr^D
1720
1721 numx1=domain_nx1*2**(refine_max_level-1)
1722 numx2=domain_nx2*2**(refine_max_level-1)
1723 numx3=domain_nx3*2**(refine_max_level-1)
1724
1725 ! parameters for creating table
1726 if (los_phi==0 .and. los_theta==90) then
1727 nxif1=domain_nx2*2**(refine_max_level-1)
1728 nxif2=domain_nx3*2**(refine_max_level-1)
1729 xifmin1=xprobmin2
1730 xifmax1=xprobmax2
1731 xifmin2=xprobmin3
1732 xifmax2=xprobmax3
1733 bnx1=block_nx2
1734 bnx2=block_nx3
1735 nbb1=domain_nx2
1736 nbb2=domain_nx3
1737 strtype1=stretch_type(2)
1738 strtype2=stretch_type(3)
1741 qs1=qstretch_baselevel(2)
1742 qs2=qstretch_baselevel(3)
1743 if (mype==0) write(*,'(a)') ' LOS vector: [-1.00 0.00 0.00]'
1744 if (mype==0) write(*,'(a)') ' xI1 vector: [ 0.00 1.00 0.00]'
1745 if (mype==0) write(*,'(a)') ' xI2 vector: [ 0.00 0.00 1.00]'
1746 else if (los_phi==90 .and. los_theta==90) then
1747 nxif1=domain_nx3*2**(refine_max_level-1)
1748 nxif2=domain_nx1*2**(refine_max_level-1)
1749 xifmin1=xprobmin3
1750 xifmax1=xprobmax3
1751 xifmin2=xprobmin1
1752 xifmax2=xprobmax1
1753 bnx1=block_nx3
1754 bnx2=block_nx1
1755 nbb1=domain_nx3
1756 nbb2=domain_nx1
1757 strtype1=stretch_type(3)
1758 strtype2=stretch_type(1)
1761 qs1=qstretch_baselevel(3)
1762 qs2=qstretch_baselevel(1)
1763 if (mype==0) write(*,'(a)') ' LOS vector: [ 0.00 -1.00 0.00]'
1764 if (mype==0) write(*,'(a)') ' xI1 vector: [-1.00 0.00 0.00]'
1765 if (mype==0) write(*,'(a)') ' xI2 vector: [ 0.00 0.00 1.00]'
1766 else
1767 nxif1=domain_nx1*2**(refine_max_level-1)
1768 nxif2=domain_nx2*2**(refine_max_level-1)
1769 xifmin1=xprobmin1
1770 xifmax1=xprobmax1
1771 xifmin2=xprobmin2
1772 xifmax2=xprobmax2
1773 bnx1=block_nx1
1774 bnx2=block_nx2
1775 nbb1=domain_nx1
1776 nbb2=domain_nx2
1777 strtype1=stretch_type(1)
1778 strtype2=stretch_type(2)
1781 qs1=qstretch_baselevel(1)
1782 qs2=qstretch_baselevel(2)
1783 if (mype==0) write(*,'(a)') ' LOS vector: [ 0.00 0.00 -1.00]'
1784 if (mype==0) write(*,'(a)') ' xI1 vector: [ 1.00 0.00 0.00]'
1785 if (mype==0) write(*,'(a)') ' xI2 vector: [ 0.00 1.00 0.00]'
1786 endif
1787 allocate(xif1(nxif1),xif2(nxif2),dxif1(nxif1),dxif2(nxif2))
1788
1789 ! initialize image coordinate
1790 select case(strtype1)
1791 case(0) ! uniform
1792 dxif1(:)=(xifmax1-xifmin1)/nxif1
1793 do ix1=1,nxif1
1794 xif1(ix1)=xifmin1+dxif1(ix1)*(ix1-half)
1795 enddo
1796 case(1) ! uni stretch
1797 qs1=qs1**(one/2**(refine_max_level-1))
1798 dxfirst1=(xifmax1-xifmin1)*(one-qs1)/(one-qs1**nxif1)
1799 dxif1(1)=dxfirst1
1800 do ix1=2,nxif1
1801 dxif1(ix1)=dxfirst1*qs1**(ix1-1)
1802 xif1(ix1)=dxif1(1)/(one-qs1)*(one-qs1**(ix1-1))+half*dxif1(ix1)
1803 enddo
1804 case(2) ! symm stretch
1805 ! base level, nbb = nstr + nuni + nstr
1806 nstr1=nstrb1*bnx1/2
1807 nuni1=nbb1-nstrb1*bnx1
1808 lenstr1=(xifmax1-xifmin1)/(2.d0+nuni1*(one-qs1)/(one-qs1**nstr1))
1809 dxfirst1=(xifmax1-xifmin1)/(dble(nuni1)+2.d0/(one-qs1)*(one-qs1**nstr1))
1810 dxmid1=dxfirst1
1811 ! refine_max level, numXI = nstr + nuni + nstr
1812 nstr1=nstr1*2**(refine_max_level-1)
1813 nuni1=nuni1*2**(refine_max_level-1)
1814 qs1=qs1**(one/2**(refine_max_level-1))
1815 dxfirst1=lenstr1*(one-qs1)/(one-qs1**nstr1)
1816 dxmid1=dxmid1/2**(refine_max_level-1)
1817 ! uniform center
1818 if(nuni1 .gt. 0) then
1819 do ix1=nstr1+1,nstr1+nuni1
1820 dxif1(ix1)=dxmid1
1821 xif1(ix1)=lenstr1+(dble(ix1)-0.5d0-nstr1)*dxif1(ix1)+xifmin1
1822 enddo
1823 endif
1824 ! left half
1825 do ix1=nstr1,1,-1
1826 dxif1(ix1)=dxfirst1*qs1**(nstr1-ix1)
1827 xif1(ix1)=xifmin1+lenstr1-dxif1(ix1)*half-dxfirst1*(one-qs1**(nstr1-ix1))/(one-qs1)
1828 enddo
1829 ! right half
1830 do ix1=nstr1+nuni1+1,nxif1
1831 dxif1(ix1)=dxfirst1*qs1**(ix1-nstr1-nuni1-1)
1832 xif1(ix1)=xifmax1-lenstr1+dxif1(ix1)*half+dxfirst1*(one-qs1**(ix1-nstr1-nuni1-1))/(one-qs1)
1833 enddo
1834 case default
1835 call mpistop("unknown stretch type")
1836 end select
1837
1838 select case(strtype2)
1839 case(0) ! uniform
1840 dxif2(:)=(xifmax2-xifmin2)/nxif2
1841 do ix2=1,nxif2
1842 xif2(ix2)=xifmin2+dxif2(ix2)*(ix2-half)
1843 enddo
1844 case(1) ! uni stretch
1845 qs2=qs2**(one/2**(refine_max_level-1))
1846 dxfirst2=(xifmax2-xifmin2)*(one-qs2)/(one-qs2**nxif2)
1847 dxif2(1)=dxfirst2
1848 do ix2=2,nxif1
1849 dxif2(ix2)=dxfirst2*qs2**(ix2-1)
1850 xif2(ix2)=dxif2(1)/(one-qs2)*(one-qs2**(ix2-1))+half*dxif2(ix2)
1851 enddo
1852 case(2) ! symm stretch
1853 ! base level, nbb = nstr + nuni + nstr
1854 nstr2=nstrb2*bnx2/2
1855 nuni2=nbb2-nstrb2*bnx2
1856 lenstr2=(xifmax2-xifmin2)/(2.d0+nuni2*(one-qs2)/(one-qs2**nstr2))
1857 dxfirst2=(xifmax2-xifmin2)/(dble(nuni2)+2.d0/(one-qs2)*(one-qs2**nstr2))
1858 dxmid2=dxfirst2
1859 ! refine_max level, numXI = nstr + nuni + nstr
1860 nstr2=nstr2*2**(refine_max_level-1)
1861 nuni2=nuni2*2**(refine_max_level-1)
1862 qs2=qs2**(one/2**(refine_max_level-1))
1863 dxfirst2=lenstr2*(one-qs2)/(one-qs2**nstr2)
1864 dxmid2=dxmid2/2**(refine_max_level-1)
1865 ! uniform center
1866 if(nuni2 .gt. 0) then
1867 do ix2=nstr2+1,nstr2+nuni2
1868 dxif2(ix2)=dxmid2
1869 xif2(ix2)=lenstr2+(dble(ix2)-0.5d0-nstr2)*dxif2(ix2)+xifmin2
1870 enddo
1871 endif
1872 ! left half
1873 do ix2=nstr2,1,-1
1874 dxif2(ix2)=dxfirst2*qs2**(nstr2-ix2)
1875 xif2(ix2)=xifmin2+lenstr2-dxif2(ix2)*half-dxfirst2*(one-qs2**(nstr2-ix2))/(one-qs2)
1876 enddo
1877 ! right half
1878 do ix2=nstr2+nuni2+1,nxif2
1879 dxif2(ix2)=dxfirst2*qs2**(ix2-nstr2-nuni2-1)
1880 xif2(ix2)=xifmax2-lenstr2+dxif2(ix2)*half+dxfirst2*(one-qs2**(ix2-nstr2-nuni2-1))/(one-qs2)
1881 enddo
1882 case default
1883 call mpistop("unknown stretch type")
1884 end select
1885
1886 ! integrate EUV flux and get cell average flux for image
1887 if (datatype=='image_euv') then
1888 if (si_unit) then
1889 unitv=unit_velocity/1.0e3 ! km/s
1890 else
1891 unitv=unit_velocity/1.0e5 ! km/s
1892 endif
1893 numwi=2
1894 allocate(wi(nxif1,nxif2,numwi))
1895 allocate(euvs(nxif1,nxif2),euv(nxif1,nxif2))
1896 allocate(dpl(nxif1,nxif2),dpls(nxif1,nxif2))
1897 euvs=0.0d0
1898 euv=0.0d0
1899 dpl=0.d0
1900 dpls=0.d0
1901 do iigrid=1,igridstail; igrid=igrids(iigrid);
1902 call integrate_euv_datresol(igrid,nxif1,nxif2,xif1,xif2,dxif1,dxif2,fl,euvs,dpls)
1903 enddo
1904 numsi=nxif1*nxif2
1905 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision, &
1906 mpi_sum,icomm,ierrmpi)
1907 call mpi_allreduce(dpls,dpl,numsi,mpi_double_precision, &
1908 mpi_sum,icomm,ierrmpi)
1909 do ix1=1,nxif1
1910 do ix2=1,nxif2
1911 if (euv(ix1,ix2)<smalldouble) euv(ix1,ix2)=zero
1912 if(euv(ix1,ix2)/=0) then
1913 dpl(ix1,ix2)=(dpl(ix1,ix2)/euv(ix1,ix2))*unitv
1914 else
1915 dpl(ix1,ix2)=0.d0
1916 endif
1917 if (abs(dpl(ix1,ix2))<smalldouble) dpl(ix1,ix2)=zero
1918 enddo
1919 enddo
1920 wi(:,:,1)=euv(:,:)
1921 wi(:,:,2)=dpl(:,:)
1922
1923 call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
1924 deallocate(wi,euv,euvs,dpl,dpls)
1925 endif
1926
1927 ! integrate SXR flux and get cell average flux for image
1928 if (datatype=='image_sxr') then
1929 if (si_unit) then
1930 arcsec=7.25d5
1931 else
1932 arcsec=7.25d7
1933 endif
1934 rhessi_rsl=2.3d0/instrument_resolution_factor
1935 numwi=1
1936 allocate(wi(nxif1,nxif2,numwi))
1937 allocate(sxrs(nxif1,nxif2),sxr(nxif1,nxif2))
1938 sxrs=0.0d0
1939 sxr=0.0d0
1940 do iigrid=1,igridstail; igrid=igrids(iigrid);
1941 call integrate_sxr_datresol(igrid,nxif1,nxif2,xif1,xif2,dxif1,dxif2,fl,sxrs)
1942 enddo
1943 numsi=nxif1*nxif2
1944 call mpi_allreduce(sxrs,sxr,numsi,mpi_double_precision, &
1945 mpi_sum,icomm,ierrmpi)
1946
1947 sxr=sxr*(rhessi_rsl*arcsec)**2 ! photons cm^-2 s^-1 pixel^-1
1948 do ix1=1,nxif1
1949 do ix2=1,nxif2
1950 if (sxr(ix1,ix2)<smalldouble) sxr(ix1,ix2)=zero
1951 enddo
1952 enddo
1953 wi(:,:,1)=sxr(:,:)
1954
1955 call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
1956 deallocate(wi,sxr,sxrs)
1957 endif
1958
1959 deallocate(xif1,xif2,dxif1,dxif2)
1960
1961 end subroutine get_image_datresol
1962
1963 subroutine integrate_sxr_datresol(igrid,nXIF1,nXIF2,xIF1,xIF2,dxIF1,dxIF2,fl,SXR)
1965
1966 integer, intent(in) :: igrid,nXIF1,nXIF2
1967 double precision, intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
1968 double precision, intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
1969 type(te_fluid), intent(in) :: fl
1970 double precision, intent(out) :: SXR(nXIF1,nXIF2)
1971
1972 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
1973 double precision :: xb^L,xd^D
1974 double precision, allocatable :: flux(:^D&)
1975 double precision, allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
1976 double precision, allocatable :: SXRg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
1977 integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
1978 double precision :: SXRt,xc^L,xg^L,r2,area_1AU
1979 integer :: ixP^L,ixP^D
1980 integer :: direction_LOS
1981
1982 if (los_phi==0 .and. los_theta==90) then
1983 direction_los=1
1984 else if (los_phi==90 .and. los_theta==90) then
1985 direction_los=2
1986 else
1987 direction_los=3
1988 endif
1989
1990 ^d&ixomin^d=ixmlo^d\
1991 ^d&ixomax^d=ixmhi^d\
1992 ^d&iximin^d=ixglo^d\
1993 ^d&iximax^d=ixghi^d\
1994 ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
1995 ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
1996
1997 allocate(flux(ixi^s))
1998 allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
1999 dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
2000 dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
2001 dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
2002 ! get local SXR flux
2003 call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr)
2004
2005 ! grid parameters
2006 levelg=ps(igrid)%level
2007 rft=2**(refine_max_level-levelg)
2008
2009 ! fine table for storing EUV flux of current grid
2010 select case(direction_los)
2011 case(1)
2012 nxg1=iximax2*rft
2013 nxg2=iximax3*rft
2014 case(2)
2015 nxg1=iximax3*rft
2016 nxg2=iximax1*rft
2017 case(3)
2018 nxg1=iximax1*rft
2019 nxg2=iximax2*rft
2020 end select
2021 allocate(sxrg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
2022 sxrg=zero
2023 xg1=zero
2024 xg2=zero
2025
2026 ! integrate for different direction
2027 select case(direction_los)
2028 case(1)
2029 do ix2=ixomin2,ixomax2
2030 ixgmin1=(ix2-1)*rft+1
2031 ixgmax1=ix2*rft
2032 do ix3=ixomin3,ixomax3
2033 ixgmin2=(ix3-1)*rft+1
2034 ixgmax2=ix3*rft
2035 sxrt=0.d0
2036 do ix1=ixomin1,ixomax1
2037 sxrt=sxrt+flux(ix^d)*dxb1(ix^d)*unit_length
2038 enddo
2039 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
2040 enddo
2041 enddo
2042 case(2)
2043 do ix3=ixomin3,ixomax3
2044 ixgmin1=(ix3-1)*rft+1
2045 ixgmax1=ix3*rft
2046 do ix1=ixomin1,ixomax1
2047 ixgmin2=(ix1-1)*rft+1
2048 ixgmax2=ix1*rft
2049 sxrt=0.d0
2050 do ix2=ixomin2,ixomax2
2051 sxrt=sxrt+flux(ix^d)*dxb2(ix^d)*unit_length
2052 enddo
2053 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
2054 enddo
2055 enddo
2056 case(3)
2057 do ix1=ixomin1,ixomax1
2058 ixgmin1=(ix1-1)*rft+1
2059 ixgmax1=ix1*rft
2060 do ix2=ixomin2,ixomax2
2061 ixgmin2=(ix2-1)*rft+1
2062 ixgmax2=ix2*rft
2063 sxrt=0.d0
2064 do ix3=ixomin3,ixomax3
2065 sxrt=sxrt+flux(ix^d)*dxb3(ix^d)*unit_length
2066 enddo
2067 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
2068 enddo
2069 enddo
2070 end select
2071
2072 area_1au=2.81d27
2073 sxrg=sxrg/area_1au
2074
2075 ! mapping grid data to global table
2076 ! index ranges in local table
2077 select case(direction_los)
2078 case(1)
2079 ixgmin1=(ixomin2-1)*rft+1
2080 ixgmax1=ixomax2*rft
2081 ixgmin2=(ixomin3-1)*rft+1
2082 ixgmax2=ixomax3*rft
2083 case(2)
2084 ixgmin1=(ixomin3-1)*rft+1
2085 ixgmax1=ixomax3*rft
2086 ixgmin2=(ixomin1-1)*rft+1
2087 ixgmax2=ixomax1*rft
2088 case(3)
2089 ixgmin1=(ixomin1-1)*rft+1
2090 ixgmax1=ixomax1*rft
2091 ixgmin2=(ixomin2-1)*rft+1
2092 ixgmax2=ixomax2*rft
2093 end select
2094 ! index ranges in global table & mapping
2095 select case(direction_los)
2096 case(1)
2097 ixpmin1=(node(pig2_,igrid)-1)*rft*block_nx2+1
2098 ixpmax1=node(pig2_,igrid)*rft*block_nx2
2099 ixpmin2=(node(pig3_,igrid)-1)*rft*block_nx3+1
2100 ixpmax2=node(pig3_,igrid)*rft*block_nx3
2101 case(2)
2102 ixpmin1=(node(pig3_,igrid)-1)*rft*block_nx3+1
2103 ixpmax1=node(pig3_,igrid)*rft*block_nx3
2104 ixpmin2=(node(pig1_,igrid)-1)*rft*block_nx1+1
2105 ixpmax2=node(pig1_,igrid)*rft*block_nx1
2106 case(3)
2107 ixpmin1=(node(pig1_,igrid)-1)*rft*block_nx1+1
2108 ixpmax1=node(pig1_,igrid)*rft*block_nx1
2109 ixpmin2=(node(pig2_,igrid)-1)*rft*block_nx2+1
2110 ixpmax2=node(pig2_,igrid)*rft*block_nx2
2111 end select
2112 xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
2113 xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
2114 dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
2115 dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
2116 sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
2117 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
2118
2119 deallocate(flux,dxb1,dxb2,dxb3,sxrg,xg1,xg2,dxg1,dxg2)
2120
2121 end subroutine integrate_sxr_datresol
2122
2123 subroutine integrate_euv_datresol(igrid,nXIF1,nXIF2,xIF1,xIF2,dxIF1,dxIF2,fl,EUV,Dpl)
2125
2126 integer, intent(in) :: igrid,nXIF1,nXIF2
2127 double precision, intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
2128 double precision, intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
2129 type(te_fluid), intent(in) :: fl
2130 double precision, intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
2131
2132 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
2133 double precision :: xb^L,xd^D
2134 double precision, allocatable :: flux(:^D&),v(:^D&),rho(:^D&)
2135 double precision, allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
2136 double precision, allocatable :: EUVg(:,:),Fvg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
2137 integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
2138 double precision :: EUVt,Fvt,xc^L,xg^L,r2
2139 integer :: ixP^L,ixP^D
2140 integer :: direction_LOS
2141
2142 if (los_phi==0 .and. los_theta==90) then
2143 direction_los=1
2144 else if (los_phi==90 .and. los_theta==90) then
2145 direction_los=2
2146 else
2147 direction_los=3
2148 endif
2149
2150 ^d&ixomin^d=ixmlo^d\
2151 ^d&ixomax^d=ixmhi^d\
2152 ^d&iximin^d=ixglo^d\
2153 ^d&iximax^d=ixghi^d\
2154 ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
2155 ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
2156
2157 allocate(flux(ixi^s),v(ixi^s),rho(ixi^s))
2158 allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
2159 dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
2160 dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
2161 dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
2162 ! get local EUV flux and velocity
2163 call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
2164 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2 ! adjust flux due to artifical change of resolution
2165 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
2166 v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
2167 deallocate(rho)
2168
2169 ! grid parameters
2170 levelg=ps(igrid)%level
2171 rft=2**(refine_max_level-levelg)
2172
2173 ! fine table for storing EUV flux of current grid
2174 select case(direction_los)
2175 case(1)
2176 nxg1=iximax2*rft
2177 nxg2=iximax3*rft
2178 case(2)
2179 nxg1=iximax3*rft
2180 nxg2=iximax1*rft
2181 case(3)
2182 nxg1=iximax1*rft
2183 nxg2=iximax2*rft
2184 end select
2185 allocate(euvg(nxg1,nxg2),fvg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
2186 euvg=zero
2187 fvg=zero
2188 xg1=zero
2189 xg2=zero
2190
2191 ! integrate for different direction
2192 select case(direction_los)
2193 case(1)
2194 do ix2=ixomin2,ixomax2
2195 ixgmin1=(ix2-1)*rft+1
2196 ixgmax1=ix2*rft
2197 do ix3=ixomin3,ixomax3
2198 ixgmin2=(ix3-1)*rft+1
2199 ixgmax2=ix3*rft
2200 euvt=0.d0
2201 fvt=0.d0
2202 do ix1=ixomin1,ixomax1
2203 euvt=euvt+flux(ix^d)*dxb1(ix^d)*unit_length
2204 fvt=fvt+flux(ix^d)*dxb1(ix^d)*unit_length*v(ix^d)
2205 enddo
2206 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
2207 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
2208 enddo
2209 enddo
2210 case(2)
2211 do ix3=ixomin3,ixomax3
2212 ixgmin1=(ix3-1)*rft+1
2213 ixgmax1=ix3*rft
2214 do ix1=ixomin1,ixomax1
2215 ixgmin2=(ix1-1)*rft+1
2216 ixgmax2=ix1*rft
2217 euvt=0.d0
2218 fvt=0.d0
2219 do ix2=ixomin2,ixomax2
2220 euvt=euvt+flux(ix^d)*dxb2(ix^d)*unit_length
2221 fvt=fvt+flux(ix^d)*dxb2(ix^d)*unit_length*v(ix^d)
2222 enddo
2223 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
2224 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
2225 enddo
2226 enddo
2227 case(3)
2228 do ix1=ixomin1,ixomax1
2229 ixgmin1=(ix1-1)*rft+1
2230 ixgmax1=ix1*rft
2231 do ix2=ixomin2,ixomax2
2232 ixgmin2=(ix2-1)*rft+1
2233 ixgmax2=ix2*rft
2234 euvt=0.d0
2235 fvt=0.d0
2236 do ix3=ixomin3,ixomax3
2237 euvt=euvt+flux(ix^d)*dxb3(ix^d)*unit_length
2238 fvt=fvt+flux(ix^d)*dxb3(ix^d)*unit_length*v(ix^d)
2239 enddo
2240 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
2241 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
2242 enddo
2243 enddo
2244 end select
2245 if (si_unit) then
2246 euvg=euvg*1.d2
2247 fvg=fvg*1.d2
2248 endif
2249
2250 ! mapping grid data to global table
2251 ! index ranges in local table
2252 select case(direction_los)
2253 case(1)
2254 ixgmin1=(ixomin2-1)*rft+1
2255 ixgmax1=ixomax2*rft
2256 ixgmin2=(ixomin3-1)*rft+1
2257 ixgmax2=ixomax3*rft
2258 case(2)
2259 ixgmin1=(ixomin3-1)*rft+1
2260 ixgmax1=ixomax3*rft
2261 ixgmin2=(ixomin1-1)*rft+1
2262 ixgmax2=ixomax1*rft
2263 case(3)
2264 ixgmin1=(ixomin1-1)*rft+1
2265 ixgmax1=ixomax1*rft
2266 ixgmin2=(ixomin2-1)*rft+1
2267 ixgmax2=ixomax2*rft
2268 end select
2269 ! index ranges in global table & mapping
2270 select case(direction_los)
2271 case(1)
2272 ixpmin1=(node(pig2_,igrid)-1)*rft*block_nx2+1
2273 ixpmax1=node(pig2_,igrid)*rft*block_nx2
2274 ixpmin2=(node(pig3_,igrid)-1)*rft*block_nx3+1
2275 ixpmax2=node(pig3_,igrid)*rft*block_nx3
2276 case(2)
2277 ixpmin1=(node(pig3_,igrid)-1)*rft*block_nx3+1
2278 ixpmax1=node(pig3_,igrid)*rft*block_nx3
2279 ixpmin2=(node(pig1_,igrid)-1)*rft*block_nx1+1
2280 ixpmax2=node(pig1_,igrid)*rft*block_nx1
2281 case(3)
2282 ixpmin1=(node(pig1_,igrid)-1)*rft*block_nx1+1
2283 ixpmax1=node(pig1_,igrid)*rft*block_nx1
2284 ixpmin2=(node(pig2_,igrid)-1)*rft*block_nx2+1
2285 ixpmax2=node(pig2_,igrid)*rft*block_nx2
2286 end select
2287 xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
2288 xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
2289 dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
2290 dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
2291 euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
2292 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
2293 dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
2294 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
2295
2296 deallocate(flux,v,dxb1,dxb2,dxb3,euvg,fvg,xg1,xg2,dxg1,dxg2)
2297
2298 end subroutine integrate_euv_datresol
2299
2300 subroutine get_image(qunit,datatype,fl)
2301 ! integrate emission flux along line of sight (LOS)
2302 ! in a 3D simulation box and get a 2D EUV image
2304 use mod_constants
2305
2306 integer, intent(in) :: qunit
2307 type(te_fluid), intent(in) :: fl
2308 character(20), intent(in) :: datatype
2309
2310 integer :: ix^D,numXI1,numXI2,numWI
2311 double precision :: xImin1,xImax1,xImin2,xImax2,xIcent1,xIcent2,dxI
2312 double precision, allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:)
2313 double precision, allocatable :: wI(:,:,:),wIs(:,:,:),EM(:,:),WLB(:,:,:)
2314 double precision :: vec_temp1(1:3),vec_temp2(1:3)
2315 double precision :: vec_z(1:3),vec_cor(1:3),xI_cor(1:2)
2316 double precision :: res,LOS_psi,r_max,r_loc
2317
2318 integer :: mass
2319 character (30) :: ion
2320 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
2321 double precision :: arcsec,RHESSI_rsl,LASCO_rsl,pixel,R_occult,smallflux
2322 integer :: iigrid,igrid,i,j,numSI,iw
2323 logical :: emit
2324
2325 if (coordinate==spherical) then
2327 else
2328 ! cartesian
2330 endif
2331
2332 ! calculate domain of the image
2333 if (coordinate==spherical) then
2334 ximin1=-abs(xprobmax1)
2335 ximin2=-abs(xprobmax1)
2336 ximax1=abs(xprobmax1)
2337 ximax2=abs(xprobmax1)
2338 else
2339 ! calculate domain of the image
2340 do ix1=1,2
2341 if (ix1==1) vec_cor(1)=xprobmin1
2342 if (ix1==2) vec_cor(1)=xprobmax1
2343 do ix2=1,2
2344 if (ix2==1) vec_cor(2)=xprobmin2
2345 if (ix2==2) vec_cor(2)=xprobmax2
2346 do ix3=1,2
2347 if (ix3==1) vec_cor(3)=xprobmin3
2348 if (ix3==2) vec_cor(3)=xprobmax3
2349 if (big_image) then
2350 r_loc=(vec_cor(1)-x_origin(1))**2
2351 r_loc=r_loc+(vec_cor(2)-x_origin(2))**2
2352 r_loc=r_loc+(vec_cor(3)-x_origin(3))**2
2353 r_loc=sqrt(r_loc)
2354 if (ix1==1 .and. ix2==1 .and. ix3==1) then
2355 r_max=r_loc
2356 else
2357 r_max=max(r_max,r_loc)
2358 endif
2359 else
2360 call get_cor_image(vec_cor,xi_cor)
2361 if (ix1==1 .and. ix2==1 .and. ix3==1) then
2362 ximin1=xi_cor(1)
2363 ximax1=xi_cor(1)
2364 ximin2=xi_cor(2)
2365 ximax2=xi_cor(2)
2366 else
2367 ximin1=min(ximin1,xi_cor(1))
2368 ximax1=max(ximax1,xi_cor(1))
2369 ximin2=min(ximin2,xi_cor(2))
2370 ximax2=max(ximax2,xi_cor(2))
2371 endif
2372 endif
2373 enddo
2374 enddo
2375 enddo
2376 if (big_image) then
2377 ximin1=-r_max
2378 ximin2=-r_max
2379 ximax1=r_max
2380 ximax2=r_max
2381 endif
2382 endif
2383 xicent1=(ximin1+ximax1)/2.d0
2384 xicent2=(ximin2+ximax2)/2.d0
2385
2386 ! tables for image
2387 if (si_unit) then
2388 arcsec=7.25d5/unit_length
2389 else
2390 arcsec=7.25d7/unit_length
2391 endif
2392 if (datatype=='image_euv') then
2393 call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
2394 dxi=spacersl*arcsec ! intrument resolution of image
2395 smallflux=smalldouble
2396 else if (datatype=='image_sxr') then
2397 rhessi_rsl=2.3d0/instrument_resolution_factor
2398 dxi=rhessi_rsl*arcsec
2399 smallflux=1.d-40
2400 else if (datatype=='image_whitelight') then
2401 if (whitelight_instrument=='LASCO/C1') then
2402 lasco_rsl=5.6d0/instrument_resolution_factor
2403 r_occult=1.1d0
2404 else if (whitelight_instrument=='LASCO/C2') then
2405 lasco_rsl=11.4d0/instrument_resolution_factor
2406 r_occult=2.d0
2407 else if (whitelight_instrument=='LASCO/C3') then
2408 lasco_rsl=56.d0/instrument_resolution_factor
2409 r_occult=3.7d0
2410 else
2411 call mpistop('Whitelight synthesis: instrument is not supported!')
2412 endif
2413 dxi=lasco_rsl*arcsec
2414 if (r_occultor>1.d0) r_occult=r_occultor
2415 r_occult=r_occult*const_rsun/unit_length
2416 smallflux=1.d-20
2417 endif
2418 numxi1=8*ceiling((ximax1-xicent1)/dxi/8.d0)
2419 ximin1=xicent1-numxi1*dxi
2420 ximax1=xicent1+numxi1*dxi
2421 numxi1=numxi1*2
2422 numxi2=8*ceiling((ximax2-xicent2)/dxi/8.d0)
2423 ximin2=xicent2-numxi2*dxi
2424 ximax2=xicent2+numxi2*dxi
2425 numxi2=numxi2*2
2426 allocate(xi1(numxi1),xi2(numxi2),dxi1(numxi1),dxi2(numxi2))
2427 do ix1=1,numxi1
2428 xi1(ix1)=ximin1+dxi*(ix1-half)
2429 dxi1(ix1)=dxi
2430 enddo
2431 do ix2=1,numxi2
2432 xi2(ix2)=ximin2+dxi*(ix2-half)
2433 dxi2(ix2)=dxi
2434 enddo
2435
2436 ! calculate emission
2437 if (datatype=='image_euv' .or. datatype=='image_sxr') then
2438 numwi=1
2439 allocate(wi(numxi1,numxi2,numwi),wis(numxi1,numxi2,numwi),em(numxi1,numxi2))
2440 wi=zero
2441 wis=zero
2442 em=zero
2443 if (coordinate==cartesian) then
2444 do iigrid=1,igridstail; igrid=igrids(iigrid);
2445 call integrate_emission_cartesian(igrid,numxi1,numxi2,xi1,xi2,dxi,fl,datatype,em)
2446 enddo
2447 else
2448 do iigrid=1,igridstail; igrid=igrids(iigrid);
2449 call integrate_emission_spherical(igrid,numxi1,numxi2,xi1,xi2,dxi,fl,datatype,em)
2450 enddo
2451 endif
2452 do ix1=1,numxi1
2453 do ix2=1,numxi2
2454 if (em(ix1,ix2)>smallflux) wis(ix1,ix2,1)=em(ix1,ix2)
2455 enddo
2456 enddo
2457 numsi=numxi1*numxi2*numwi
2458 call mpi_allreduce(wis,wi,numsi,mpi_double_precision,mpi_sum,icomm,ierrmpi)
2459 if (activate_unit_arcsec) then
2460 xi1=xi1/arcsec
2461 dxi1=dxi1/arcsec
2462 xi2=xi2/arcsec
2463 dxi2=dxi2/arcsec
2464 endif
2465 call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
2466 deallocate(wi,wis,em)
2467 else if (datatype=='image_whitelight') then
2468 numwi=2
2469 allocate(wi(numxi1,numxi2,numwi),wis(numxi1,numxi2,numwi),wlb(numxi1,numxi2,numwi))
2470 wi=zero
2471 wis=zero
2472 wlb=zero
2473 if (coordinate==spherical) then
2474 do iigrid=1,igridstail; igrid=igrids(iigrid);
2475 call integrate_whitelight_spherical(igrid,numxi1,numxi2,numwi,xi1,xi2,dxi,fl,datatype,wlb)
2476 enddo
2477 endif
2478 do ix1=1,numxi1
2479 do ix2=1,numxi2
2480 if (wlb(ix1,ix2,1)>smallflux) then
2481 wis(ix1,ix2,1)=wlb(ix1,ix2,1)
2482 wis(ix1,ix2,2)=wlb(ix1,ix2,2)
2483 endif
2484 enddo
2485 enddo
2486 numsi=numxi1*numxi2*numwi
2487 call mpi_allreduce(wis,wi,numsi,mpi_double_precision,mpi_sum,icomm,ierrmpi)
2488 if (activate_unit_arcsec) then
2489 xi1=xi1/arcsec
2490 dxi1=dxi1/arcsec
2491 xi2=xi2/arcsec
2492 dxi2=dxi2/arcsec
2493 endif
2494 call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
2495 deallocate(wi,wis,wlb)
2496 endif
2497
2498 deallocate(xi1,xi2,dxi1,dxi2)
2499
2500 end subroutine get_image
2501
2502 subroutine integrate_emission_cartesian(igrid,numXI1,numXI2,xI1,xI2,dxI,fl,datatype,EM)
2503 integer, intent(in) :: igrid,numXI1,numXI2
2504 double precision, intent(in) :: xI1(numXI1),xI2(numXI2)
2505 double precision, intent(in) :: dxI
2506 type(te_fluid), intent(in) :: fl
2507 character(20), intent(in) :: datatype
2508 double precision, intent(inout) :: EM(numXI1,numXI2)
2509
2510 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
2511 double precision :: xb^L,xd^D
2512 double precision, allocatable :: flux(:^D&)
2513 double precision :: res
2514 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
2515 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC
2516 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2)
2517
2518 integer :: mass
2519 double precision :: logTe
2520 character (30) :: ion
2521 double precision :: lineCent
2522 double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
2523 double precision :: arcsec,pixel,RHESSI_rsl,area_1AU
2524 double precision :: aa,bb
2525
2526 ^d&ixomin^d=ixmlo^d\
2527 ^d&ixomax^d=ixmhi^d\
2528 ^d&iximin^d=ixglo^d\
2529 ^d&iximax^d=ixghi^d\
2530 ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
2531 ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
2532
2533 if (si_unit) then
2534 arcsec=7.25d5/unit_length
2535 else
2536 arcsec=7.25d7/unit_length
2537 endif
2538
2539 allocate(flux(ixi^s))
2540 if (datatype=='image_euv') then
2541 ! get local EUV flux and velocity
2542 call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
2543 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2 ! adjust flux due to artifical change of resolution
2544 call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
2545 pixel=spacersl*arcsec
2546 sigma0=sigma_psf*pixel
2547 else if (datatype=='image_sxr') then
2548 ! get local SXR flux photons cm^-3 s^-1 (cgs) or photons m^-3 s^-1 (SI)
2549 call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr)
2550 rhessi_rsl=2.3d0/instrument_resolution_factor
2551 sigma_psf=1.d0
2552 pixel=rhessi_rsl*arcsec
2553 sigma0=sigma_psf*pixel
2554 area_1au=2.81d27
2555 endif
2556
2557 ! integrate emission
2558 {do ix^d=ixomin^d,ixomax^d\}
2559 ^d&nsubc^d=1;
2560 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(vec_xi1(^d))/(dxi/2.d0)));
2561 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(vec_xi2(^d))/(dxi/2.d0)));
2562 ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
2563 if (datatype=='image_euv') then
2564 if (si_unit) then
2565 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length*1.d2/dxi/dxi ! DN s^-1
2566 else
2567 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length/dxi/dxi ! DN s^-1
2568 endif
2569 else if (datatype=='image_sxr') then
2570 ! sub-cell SXR flux at 1 AU [photons s^-1 cm^-2]
2571 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length**3/area_1au
2572 endif
2573 if (fluxsubc>smalldouble) then
2574 ! dividing a cell to several parts to get more accurate integrating values
2575 {do isubc^d=1,nsubc^d\}
2576 ^d&xsubc(^d)=ps(igrid)%x(ix^dd,^d)-half*ps(igrid)%dx(ix^dd,^d)+(isubc^d-half)*dxsubc^d;
2577 ! mapping the 3D coordinate to location at the image
2578 call get_cor_image(xsubc,xcent)
2579 ! distribution at nearby pixels
2580 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
2581 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
2582 ixpmin1=max(1,ixp1-3)
2583 ixpmax1=min(ixp1+3,numxi1)
2584 ixpmin2=max(1,ixp2-3)
2585 ixpmax2=min(ixp2+3,numxi2)
2586 do ixp1=ixpmin1,ixpmax1
2587 do ixp2=ixpmin2,ixpmax2
2588 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
2589 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
2590 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
2591 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
2592 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
2593 em(ixp1,ixp2)=em(ixp1,ixp2)+fluxsubc*factor
2594 enddo !ixP2
2595 enddo !ixP1
2596 {enddo\} !iSubC
2597 endif
2598 {enddo\} !ix
2599
2600 deallocate(flux)
2601 end subroutine integrate_emission_cartesian
2602
2603 subroutine integrate_emission_spherical(igrid,numXI1,numXI2,xI1,xI2,dxI,fl,datatype,EM)
2604 integer, intent(in) :: igrid,numXI1,numXI2
2605 double precision, intent(in) :: xI1(numXI1),xI2(numXI2)
2606 double precision, intent(in) :: dxI
2607 type(te_fluid), intent(in) :: fl
2608 character(20), intent(in) :: datatype
2609 double precision, intent(inout) :: EM(numXI1,numXI2)
2610
2611 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
2612 double precision, allocatable :: flux(:^D&),Ne(:^D&)
2613 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
2614 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC,RsubC
2615 double precision :: TBsubC,PBsubC
2616 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2),xSubC_car(1:3)
2617 double precision :: R_thick,dotp,dvolume,R_occult,Rc
2618 double precision :: dxl(1:3),x_sph(1:3),dx_sph(1:3)
2619 double precision :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
2620 logical :: sun_back_side,emit
2621
2622 integer :: mass
2623 double precision :: logTe
2624 character (30) :: ion
2625 double precision :: lineCent
2626 double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
2627 double precision :: RHESSI_rsl,area_1AU,arcsec,pixel
2628
2629 ^d&ixomin^d=ixmlo^d;
2630 ^d&ixomax^d=ixmhi^d;
2631 ^d&iximin^d=ixglo^d;
2632 ^d&iximax^d=ixghi^d;
2633
2634 if (si_unit) then
2635 arcsec=7.25d5/unit_length
2636 else
2637 arcsec=7.25d7/unit_length
2638 endif
2639
2640 allocate(flux(ixi^s))
2641 if (datatype=='image_euv') then
2642 ! get local EUV flux and velocity
2643 call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
2644 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2 ! adjust flux due to artifical change of resolution
2645 call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
2646 pixel=spacersl*arcsec
2647 sigma0=sigma_psf*pixel
2648 else if (datatype=='image_sxr') then
2649 ! get local SXR flux photons cm^-3 s^-1 (cgs) or photons m^-3 s^-1 (SI)
2650 call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr)
2651 rhessi_rsl=2.3d0/instrument_resolution_factor
2652 sigma_psf=1.d0
2653 pixel=rhessi_rsl*arcsec
2654 sigma0=sigma_psf*pixel
2655 area_1au=2.81d27
2656 endif
2657
2658 ! integrate emission
2659 r_thick=r_opt_thick*const_rsun/unit_length
2660 {do ix^d=ixomin^d,ixomax^d\}
2661 x_sph(1:3)=ps(igrid)%x(ix^d,1:3)
2662 dx_sph(1:3)=ps(igrid)%dx(ix^d,1:3)
2663 dxl(1)=dx_sph(1) ! cell size in length
2664 dxl(2)=x_sph(1)*dx_sph(2) ! cell size in length
2665 dxl(3)=x_sph(1)*dsin(x_sph(2))*dx_sph(3) ! cell size in length
2666 ! dividing a cell to several sub-cells to get more accurate integrating values
2667 ^d&nsubc^d=1;
2668 call get_unit_vector_spherical(x_sph,unitv_r,unitv_theta,unitv_phi)
2669 call dot_product_loc(unitv_r,vec_xi1,dotp)
2670 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
2671 call dot_product_loc(unitv_r,vec_xi2,dotp)
2672 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
2673 call dot_product_loc(unitv_theta,vec_xi1,dotp)
2674 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
2675 call dot_product_loc(unitv_theta,vec_xi2,dotp)
2676 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
2677 call dot_product_loc(unitv_phi,vec_xi1,dotp)
2678 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
2679 call dot_product_loc(unitv_phi,vec_xi2,dotp)
2680 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
2681
2682 ! integrate sub-cells
2683 do isubc1=1,nsubc1
2684 ! sub-cell center coordinate in spherical
2685 xsubc(1)=x_sph(1)-half*dx_sph(1)+(isubc1-half)*dx_sph(1)/nsubc1
2686 rsubc=xsubc(1)
2687 dxsubc1=dx_sph(1)/nsubc1 ! sub-cell size in length
2688 do isubc2=1,nsubc2
2689 ! sub-cell center coordinate in spherical
2690 xsubc(2)=x_sph(2)-half*dx_sph(2)+(isubc2-half)*dx_sph(2)/nsubc2
2691 dxsubc2=xsubc(1)*dx_sph(2)/nsubc2 ! sub-cell size in length
2692 dxsubc3=xsubc(1)*dsin(xsubc(2))*dx_sph(3)/nsubc3 ! sub-cell size in length
2693 dvolume=dxsubc1*dxsubc2*dxsubc3
2694 if (datatype=='image_euv') then
2695 if (si_unit) then
2696 fluxsubc=flux(ix^d)*dvolume*unit_length*1.d2/dxi/dxi ! DN s^-1
2697 else
2698 fluxsubc=flux(ix^d)*dvolume*unit_length/dxi/dxi ! DN s^-1
2699 endif
2700 else if (datatype=='image_sxr') then
2701 ! sub-cell SXR flux at 1 AU [photons s^-1 cm^-2]
2702 fluxsubc=flux(ix^d)*dvolume*unit_length**3/area_1au
2703 endif
2704 ! enter integration if flux large enough
2705 if (fluxsubc>smalldouble) then
2706 do isubc3=1,nsubc3
2707 ! sub-cell center coordinate in spherical
2708 xsubc(3)=x_sph(3)-half*dx_sph(3)+(isubc3-half)*dx_sph(3)/nsubc3
2709 call get_cor_image_spherical(xsubc,xcent)
2710 rc=dsqrt(xcent(1)**2+xcent(2)**2) ! distance to sun center (on the image plane)
2711 !
2712 ! whether the local emitted photons can arrive the telescope
2713 call spherical_to_cartesian(xsubc,xsubc_car)
2714 call dot_product_loc(vec_los,xsubc_car,dotp)
2715 sun_back_side=.true.
2716 if (dotp<0.d0) sun_back_side=.false.
2717 ! whether the local emission can reach the telescope
2718 if (sun_back_side) then
2719 emit=.false.
2720 if (rc>r_thick) emit=.true.
2721 else
2722 emit=.true.
2723 if (xsubc(1)<=r_thick) emit=.false.
2724 endif
2725 !
2726 if (emit) then
2727 ! mapping the 3D coordinate to location at the image
2728 ! distribution at nearby pixels
2729 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
2730 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
2731 ixpmin1=max(1,ixp1-3)
2732 ixpmax1=min(ixp1+3,numxi1)
2733 ixpmin2=max(1,ixp2-3)
2734 ixpmax2=min(ixp2+3,numxi2)
2735 do ixp1=ixpmin1,ixpmax1
2736 do ixp2=ixpmin2,ixpmax2
2737 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
2738 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
2739 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
2740 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
2741 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
2742 em(ixp1,ixp2)=em(ixp1,ixp2)+fluxsubc*factor
2743 enddo !ixP2
2744 enddo !ixP1
2745 endif !emit
2746 enddo !iSubC3
2747 endif !smallflux
2748 enddo !iSubC2
2749 enddo !iSubC1
2750 {enddo\} !ix
2751
2752 deallocate(flux)
2753
2754 end subroutine integrate_emission_spherical
2755
2756 subroutine integrate_whitelight_spherical(igrid,numXI1,numXI2,numWI,xI1,xI2,dxI,fl,datatype,WLB)
2757 integer, intent(in) :: igrid,numXI1,numXI2,numWI
2758 double precision, intent(in) :: xI1(numXI1),xI2(numXI2)
2759 double precision, intent(in) :: dxI
2760 type(te_fluid), intent(in) :: fl
2761 character(20), intent(in) :: datatype
2762 double precision, intent(inout) :: WLB(numXI1,numXI2,numWI)
2763
2764 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
2765 double precision, allocatable :: flux(:^D&),Ne(:^D&)
2766 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
2767 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC,RsubC
2768 double precision :: sigma_PSF,sigma0,arcsec,pixel,LASCO_rsl
2769 double precision :: A,B,C,D,Rc,Ne0,TBsubC,PBsubC,factor
2770 double precision :: R_thick,dotp,dvolume,R_occult
2771 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2),xSubC_car(1:3)
2772 double precision :: dxl(1:3),x_sph(1:3),dx_sph(1:3)
2773 double precision :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
2774 logical :: emit
2775
2776 ^d&ixomin^d=ixmlo^d;
2777 ^d&ixomax^d=ixmhi^d;
2778 ^d&iximin^d=ixglo^d;
2779 ^d&iximax^d=ixghi^d;
2780
2781 if (si_unit) then
2782 arcsec=7.25d5/unit_length
2783 else
2784 arcsec=7.25d7/unit_length
2785 endif
2786
2787 allocate(ne(ixi^s))
2788 if (whitelight_instrument=='LASCO/C1') then
2789 lasco_rsl=5.6d0/instrument_resolution_factor
2790 r_occult=1.1d0
2791 else if (whitelight_instrument=='LASCO/C2') then
2792 lasco_rsl=11.4d0/instrument_resolution_factor
2793 r_occult=2.d0
2794 else if (whitelight_instrument=='LASCO/C3') then
2795 lasco_rsl=56.d0/instrument_resolution_factor
2796 r_occult=3.7d0
2797 endif
2798 if (r_occultor>1.d0) r_occult=r_occultor
2799 r_occult=r_occult*const_rsun/unit_length
2800 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,ne)
2801 sigma_psf=1.d0
2802 pixel=lasco_rsl*arcsec
2803 sigma0=sigma_psf*pixel
2804
2805 ! integrate emission
2806 r_thick=r_opt_thick*const_rsun/unit_length
2807 {do ix^d=ixomin^d,ixomax^d\}
2808 x_sph(1:3)=ps(igrid)%x(ix^d,1:3)
2809 dx_sph(1:3)=ps(igrid)%dx(ix^d,1:3)
2810 dxl(1)=dx_sph(1) ! cell size in length
2811 dxl(2)=x_sph(1)*dx_sph(2) ! cell size in length
2812 dxl(3)=x_sph(1)*dsin(x_sph(2))*dx_sph(3) ! cell size in length
2813 ne0=ne(ix^d)*unit_numberdensity
2814 ! dividing a cell to several sub-cells to get more accurate integrating values
2815 ^d&nsubc^d=1;
2816 call get_unit_vector_spherical(x_sph,unitv_r,unitv_theta,unitv_phi)
2817 call dot_product_loc(unitv_r,vec_xi1,dotp)
2818 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
2819 call dot_product_loc(unitv_r,vec_xi2,dotp)
2820 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
2821 call dot_product_loc(unitv_theta,vec_xi1,dotp)
2822 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
2823 call dot_product_loc(unitv_theta,vec_xi2,dotp)
2824 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
2825 call dot_product_loc(unitv_phi,vec_xi1,dotp)
2826 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
2827 call dot_product_loc(unitv_phi,vec_xi2,dotp)
2828 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
2829
2830 ! integrate sub-cells
2831 do isubc1=1,nsubc1
2832 ! sub-cell center coordinate in spherical
2833 xsubc(1)=x_sph(1)-half*dx_sph(1)+(isubc1-half)*dx_sph(1)/nsubc1
2834 rsubc=xsubc(1)
2835 dxsubc1=dx_sph(1)/nsubc1 ! sub-cell size in length
2836 call get_thomson_parameters(rsubc,a,b,c,d)
2837 do isubc2=1,nsubc2
2838 ! sub-cell center coordinate in spherical
2839 xsubc(2)=x_sph(2)-half*dx_sph(2)+(isubc2-half)*dx_sph(2)/nsubc2
2840 dxsubc2=xsubc(1)*dx_sph(2)/nsubc2 ! sub-cell size in length
2841 dxsubc3=xsubc(1)*dsin(xsubc(2))*dx_sph(3)/nsubc3 ! sub-cell size in length
2842 dvolume=dxsubc1*dxsubc2*dxsubc3
2843 do isubc3=1,nsubc3
2844 ! sub-cell center coordinate in spherical
2845 xsubc(3)=x_sph(3)-half*dx_sph(3)+(isubc3-half)*dx_sph(3)/nsubc3
2846 call get_cor_image_spherical(xsubc,xcent)
2847 rc=dsqrt(xcent(1)**2+xcent(2)**2) ! distance to sun center (on the image plane)
2848 ! whether the local emitted photons can arrive the telescope
2849 emit=.false.
2850 if (rc>r_occult) then
2851 emit=.true.
2852 ! scaterring flux from cm^-3 of plasma
2853 call get_whitelight_thomson(rsubc,rc,ne0,a,b,c,d,tbsubc,pbsubc)
2854 tbsubc=tbsubc*dvolume*unit_length/dxi/dxi
2855 pbsubc=pbsubc*dvolume*unit_length/dxi/dxi
2856 if (tbsubc<1.d-20) emit=.false.
2857 endif
2858 if (emit) then
2859 ! mapping the 3D coordinate to location at the image
2860 ! distribution at nearby pixels
2861 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
2862 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
2863 ixpmin1=max(1,ixp1-3)
2864 ixpmax1=min(ixp1+3,numxi1)
2865 ixpmin2=max(1,ixp2-3)
2866 ixpmax2=min(ixp2+3,numxi2)
2867 do ixp1=ixpmin1,ixpmax1
2868 do ixp2=ixpmin2,ixpmax2
2869 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
2870 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
2871 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
2872 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
2873 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
2874 wlb(ixp1,ixp2,1)=wlb(ixp1,ixp2,1)+tbsubc*factor
2875 wlb(ixp1,ixp2,2)=wlb(ixp1,ixp2,2)+pbsubc*factor
2876 enddo !ixP2
2877 enddo !ixP1
2878 endif
2879 enddo !iSubC3
2880 enddo !iSubC2
2881 enddo !iSubC1
2882 {enddo\} !ix
2883
2884 deallocate(ne)
2885
2886 end subroutine integrate_whitelight_spherical
2887
2888 subroutine get_thomson_parameters(Rl,A,B,C,D)
2889 ! parameters given in Billings 1968
2890 use mod_constants
2891 double precision, intent(in) :: Rl
2892 double precision, intent(inout) :: A,B,C,D
2893
2894 double precision :: sinO,cosO,sinO2,cosO2,tmp
2895
2896 sino=const_rsun/(rl*unit_length)
2897 sino2=sino**2
2898 coso2=1.d0-sino2
2899 coso=abs(dsqrt(coso2))
2900 tmp=log((1.d0+sino)/coso)
2901 a=coso*sino2
2902 b=-(1.d0-3.d0*sino2-(coso2/sino)*(1.d0+3.d0*sino2)*tmp)/8.d0
2903 c=4.d0/3.d0-coso-coso*coso2/3.d0
2904 d=(5.d0+sino2-(coso2/sino)*(5.d0-sino2)*tmp)/8.d0
2905
2906 end subroutine get_thomson_parameters
2907
2908 subroutine get_whitelight_thomson(Rl,Rin,Ne,A,B,C,D,fluxTB,fluxPB)
2909 ! use the method in SSW/eltheory
2910 double precision, intent(in) :: Rl,Rin,Ne,A,B,C,D
2911 double precision, intent(inout) :: fluxTB,fluxPB
2912
2913 double precision :: const,u,Bt,Br,PB,TB,sinchi2
2914
2915 u=0.63d0
2916 const=1.24878d-25/(1.d0-u/3.d0)
2917 sinchi2=(rin/rl)**2
2918 bt=const*(c+u*(d-c))
2919 pb=const*sinchi2*((a+u*(b-a)))
2920 br=bt-pb
2921 tb=bt+br
2922 fluxtb=tb*ne
2923 fluxpb=pb*ne
2924
2925 end subroutine get_whitelight_thomson
2926
2927 subroutine get_unit_vector_spherical(x_sph,unitv_r,unitv_theta,unitv_phi)
2928 double precision, intent(in) :: x_sph(1:3)
2929 double precision, intent(inout) :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
2930
2931 unitv_r(1)=dsin(x_sph(2))*dcos(x_sph(3))
2932 unitv_r(2)=dsin(x_sph(2))*dsin(x_sph(3))
2933 unitv_r(3)=dcos(x_sph(2))
2934 unitv_theta(1)=dcos(x_sph(2))*dcos(x_sph(3))
2935 unitv_theta(2)=dcos(x_sph(2))*dsin(x_sph(3))
2936 unitv_theta(3)=-dsin(x_sph(2))
2937 unitv_phi(1)=-dsin(x_sph(3))
2938 unitv_phi(2)=dcos(x_sph(3))
2939 unitv_phi(3)=zero
2940
2941 end subroutine get_unit_vector_spherical
2942
2943 subroutine output_data(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,datatype)
2944 ! change the format of data and write data
2946
2947 integer, intent(in) :: qunit,nXO1,nXO2,nWO
2948 double precision, intent(in) :: dxO1(nxO1),dxO2(nxO2)
2949 double precision, intent(in) :: xO1(nXO1),xO2(nxO2)
2950 double precision, intent(inout) :: wO(nXO1,nXO2,nWO)
2951 character(20), intent(in) :: datatype
2952
2953 integer :: nPiece,nP1,nP2,nC1,nC2,nWC
2954 integer :: piece_nmax1,piece_nmax2,ix1,ix2,j,ipc,ixc1,ixc2
2955 double precision, allocatable :: xC(:,:,:,:),wC(:,:,:,:),dxC(:,:,:,:)
2956
2957 ! clean small values
2958 do ix1=1,nxo1
2959 do ix2=1,nxo2
2960 do j=1,nwo
2961 if (abs(wo(ix1,ix2,j))<smalldouble) wo(ix1,ix2,j)=zero
2962 enddo
2963 enddo
2964 enddo
2965
2966 ! how many cells in each grid
2967 if (dat_resolution) then
2968 if (datatype=='image_euv' .or. datatype=='image_sxr') then
2969 if (los_phi==0 .and. los_theta==90) then
2970 piece_nmax1=block_nx2
2971 piece_nmax2=block_nx3
2972 else if (los_phi==90 .and. los_theta==90) then
2973 piece_nmax1=block_nx3
2974 piece_nmax2=block_nx1
2975 else
2976 piece_nmax1=block_nx1
2977 piece_nmax2=block_nx2
2978 endif
2979 else if (datatype=='spectrum_euv') then
2980 piece_nmax1=16
2981 if (direction_slit==1) then
2982 piece_nmax2=block_nx1
2983 else if (direction_slit==2) then
2984 piece_nmax2=block_nx2
2985 else
2986 piece_nmax2=block_nx3
2987 endif
2988 endif
2989 else
2990 piece_nmax1=20
2991 piece_nmax2=20
2992 endif
2993 loopn1: do j=piece_nmax1,1,-1
2994 if(mod(nxo1,j)==0) then
2995 nc1=j
2996 exit loopn1
2997 endif
2998 enddo loopn1
2999 loopn2: do j=piece_nmax2,1,-1
3000 if(mod(nxo2,j)==0) then
3001 nc2=j
3002 exit loopn2
3003 endif
3004 enddo loopn2
3005 ! how many grids
3006 np1=nxo1/nc1
3007 np2=nxo2/nc2
3008 npiece=np1*np2
3009 nwc=nwo
3010
3011 ! output images
3012 select case(convert_type)
3013 case('EIvtuCCmpi','ESvtuCCmpi','SIvtuCCmpi','WIvtuCCmpi')
3014 ! put data into grids
3015 allocate(xc(npiece,nc1,nc2,2))
3016 allocate(dxc(npiece,nc1,nc2,2))
3017 allocate(wc(npiece,nc1,nc2,nwo))
3018 do ipc=1,npiece
3019 do ixc1=1,nc1
3020 do ixc2=1,nc2
3021 ix1=mod(ipc-1,np1)*nc1+ixc1
3022 ix2=floor(1.0*(ipc-1)/np1)*nc2+ixc2
3023 xc(ipc,ixc1,ixc2,1)=xo1(ix1)
3024 xc(ipc,ixc1,ixc2,2)=xo2(ix2)
3025 dxc(ipc,ixc1,ixc2,1)=dxo1(ix1)
3026 dxc(ipc,ixc1,ixc2,2)=dxo2(ix2)
3027 do j=1,nwc
3028 wc(ipc,ixc1,ixc2,j)=wo(ix1,ix2,j)
3029 enddo
3030 enddo
3031 enddo
3032 enddo
3033 ! write data into vtu file
3034 call write_image_vtucc(qunit,xc,wc,dxc,npiece,nc1,nc2,nwc,datatype)
3035 deallocate(xc,dxc,wc)
3036 case('EIvtiCCmpi','ESvtiCCmpi','SIvtiCCmpi','WIvtiCCmpi')
3037 if (sum(stretch_type(:))>0 .and. dat_resolution) then
3038 call mpistop("Error in synthesize emission: vti is not supported for dat resolution")
3039 else
3040 call write_image_vticc(qunit,xo1,xo2,dxo1,dxo2,wo,nxo1,nxo2,nwo,nc1,nc2)
3041 endif
3042 case default
3043 call mpistop("Error in synthesize emission: Unknown convert_type")
3044 end select
3045
3046 end subroutine output_data
3047 }
3048
3049 subroutine write_image_vticc(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,nC1,nC2)
3050 ! write image data to vti
3052
3053 integer, intent(in) :: qunit,nXO1,nXO2,nWO,nC1,nC2
3054 double precision, intent(in) :: xO1(nXO1),xO2(nxO2)
3055 double precision, intent(in) :: dxO1(nxO1),dxO2(nxO2)
3056 double precision, intent(in) :: wO(nXO1,nXO2,nWO)
3057
3058 double precision :: origin(1:3), spacing(1:3)
3059 integer :: wholeExtent(1:6),extent(1:6)
3060 integer :: nP1,nP2,iP1,iP2,iw
3061 integer :: ixC1,ixC2,ixCmin1,ixCmax1,ixCmin2,ixCmax2
3062
3063 integer :: filenr
3064 logical :: fileopen
3065 character (70) :: subname,wname,vname,nameL,nameS
3066 character (len=std_len) :: filename
3067 integer :: mass
3068 double precision :: logTe
3069 character (30) :: ion
3070 double precision :: line_center
3071 double precision :: spatial_rsl,spectral_rsl,sigma_PSF,wslit
3072
3073
3074 origin(1)=xo1(1)-0.5d0*dxo1(1)
3075 origin(2)=xo2(1)-0.5d0*dxo2(1)
3076 origin(3)=zero
3077 spacing(1)=dxo1(1)
3078 spacing(2)=dxo2(1)
3079 spacing(3)=zero
3080 wholeextent=0
3081 wholeextent(2)=nxo1
3082 wholeextent(4)=nxo2
3083 np1=nxo1/nc1
3084 np2=nxo2/nc2
3085
3086 ! get information of emission line
3087 if (convert_type=='EIvtiCCmpi') then
3088 call get_line_info(wavelength,ion,mass,logte,line_center,spatial_rsl,spectral_rsl,sigma_psf,wslit)
3089 else if (convert_type=='ESvtiCCmpi') then
3090 call get_line_info(spectrum_wl,ion,mass,logte,line_center,spatial_rsl,spectral_rsl,sigma_psf,wslit)
3091 endif
3092
3093 if (mype==0) then
3094 inquire(qunit,opened=fileopen)
3095 if(.not.fileopen)then
3096 ! generate filename
3097 filenr=snapshotini
3098 if (autoconvert) filenr=snapshotnext
3099 if (convert_type=='EIvtiCCmpi') then
3100 write(filename,'(a,i4.4,a)') trim(filename_euv),filenr,".vti"
3101 else if (convert_type=='SIvtiCCmpi') then
3102 write(filename,'(a,i4.4,a)') trim(filename_sxr),filenr,".vti"
3103 else if (convert_type=='WIvtiCCmpi') then
3104 write(filename,'(a,i4.4,a)') trim(filename_whitelight),filenr,".vti"
3105 else if (convert_type=='ESvtiCCmpi') then
3106 write(filename,'(a,i4.4,a)') trim(filename_spectrum),filenr,".vti"
3107 endif
3108 open(qunit,file=filename,status='unknown',form='formatted')
3109 endif
3110
3111 ! generate xml header
3112 write(qunit,'(a)')'<?xml version="1.0"?>'
3113 write(qunit,'(a)',advance='no') '<VTKFile type="ImageData"'
3114 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
3115 write(qunit,'(a,3(1pe14.6),a,6(i10),a,3(1pe14.6),a)')' <ImageData Origin="',&
3116 origin,'" WholeExtent="',wholeextent,'" Spacing="',spacing,'">'
3117 ! file info
3118 write(qunit,'(a)')'<FieldData>'
3119 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
3120 'NumberOfTuples="1" format="ascii">'
3121 write(qunit,*) real(global_time*time_convert_factor)
3122 write(qunit,'(a)')'</DataArray>'
3123 if (convert_type=='EIvtiCCmpi' .or. convert_type=='ESvtiCCmpi') then
3124 write(qunit,'(2a)')'<DataArray type="Float32" Name="logT" ',&
3125 'NumberOfTuples="1" format="ascii">'
3126 write(qunit,*) real(logte)
3127 write(qunit,'(a)')'</DataArray>'
3128 endif
3129 write(qunit,'(a)')'</FieldData>'
3130 ! pixel/cell data
3131 do ip1=1,np1
3132 do ip2=1,np2
3133 extent=0
3134 extent(1)=(ip1-1)*nc1
3135 extent(2)=ip1*nc1
3136 extent(3)=(ip2-1)*nc2
3137 extent(4)=ip2*nc2
3138 ixcmin1=extent(1)+1
3139 ixcmax1=extent(2)
3140 ixcmin2=extent(3)+1
3141 ixcmax2=extent(4)
3142 write(qunit,'(a,6(i10),a)') &
3143 '<Piece Extent="',extent,'">'
3144 write(qunit,'(a)')'<CellData>'
3145 do iw=1,nwo
3146 ! variable name
3147 if (convert_type=='EIvtiCCmpi') then
3148 if (wavelength<100) then
3149 write(vname,'(a,i2)') "AIA",wavelength
3150 else if (wavelength<1000) then
3151 write(vname,'(a,i3)') "AIA",wavelength
3152 else
3153 write(vname,'(a,i4)') "IRIS",wavelength
3154 endif
3155 if (iw==2 .and. dat_resolution) vname='Doppler_velocity'
3156 else if (convert_type=='SIvtiCCmpi') then
3157 if (emin_sxr<10 .and. emax_sxr<10) then
3158 write(vname,'(a,i1,a,i1,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
3159 else if (emin_sxr<10 .and. emax_sxr>=10) then
3160 write(vname,'(a,i1,a,i2,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
3161 else
3162 write(vname,'(a,i2,a,i2,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
3163 endif
3164 else if (convert_type=='WIvtiCCmpi') then
3165 if (iw==1) write(vname,'(a)')'B'
3166 if (iw==2) write(vname,'(a)')'pB'
3167 else if (convert_type=='ESvtiCCmpi') then
3168 if (spectrum_wl==1354) then
3169 write(vname,'(a,i4)') "SG",spectrum_wl
3170 else
3171 write(vname,'(a,i3)') "EIS",spectrum_wl
3172 endif
3173 endif
3174 write(qunit,'(a,a,a)')&
3175 '<DataArray type="Float64" Name="',trim(vname),'" format="ascii">'
3176 write(qunit,'(200(1pe14.6))') ((wo(ixc1,ixc2,iw),ixc1=ixcmin1,ixcmax1),ixc2=ixcmin2,ixcmax2)
3177 write(qunit,'(a)')'</DataArray>'
3178 enddo
3179 write(qunit,'(a)')'</CellData>'
3180 write(qunit,'(a)')'</Piece>'
3181 enddo
3182 enddo
3183 ! end
3184 write(qunit,'(a)')'</ImageData>'
3185 write(qunit,'(a)')'</VTKFile>'
3186 close(qunit)
3187 endif
3188
3189 end subroutine write_image_vticc
3190
3191 subroutine write_image_vtucc(qunit,xC,wC,dxC,nPiece,nC1,nC2,nWC,datatype)
3192 ! write image data to vtu
3194
3195 integer, intent(in) :: qunit
3196 integer, intent(in) :: nPiece,nC1,nC2,nWC
3197 double precision, intent(in) :: xC(nPiece,nC1,nC2,2),dxC(nPiece,nc1,nc2,2)
3198 double precision, intent(in) :: wC(nPiece,nC1,nC2,nWC)
3199 character(20), intent(in) :: datatype
3200
3201 integer :: nP1,nP2
3202 double precision :: xP(nPiece,nC1+1,nC2+1,2)
3203 integer :: filenr
3204 logical :: fileopen
3205 character (70) :: subname,wname,vname,nameL,nameS
3206 character (len=std_len) :: filename
3207 integer :: ixC1,ixC2,ixP,ix1,ix2,j
3208 integer :: nc,np,icel,VTK_type
3209 integer :: mass
3210 double precision :: logTe
3211 character (30) :: ion
3212 double precision :: line_center
3213 double precision :: spatial_rsl,spectral_rsl,sigma_PSF,wslit
3214
3215 np1=nc1+1
3216 np2=nc2+1
3217 np=np1*np2
3218 nc=nc1*nc2
3219 ! cell corner location
3220 do ixp=1,npiece
3221 do ix1=1,np1
3222 do ix2=1,np2
3223 if (ix1<np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1,1,1)-0.5d0*dxc(ixp,ix1,1,1)
3224 if (ix1==np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1-1,1,1)+0.5d0*dxc(ixp,ix1-1,1,1)
3225 if (ix2<np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2,2)-0.5d0*dxc(ixp,1,ix2,2)
3226 if (ix2==np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2-1,2)+0.5d0*dxc(ixp,1,ix2-1,2)
3227 enddo
3228 enddo
3229 enddo
3230 ! get information of emission line
3231 if (datatype=='image_euv') then
3232 call get_line_info(wavelength,ion,mass,logte,line_center,spatial_rsl,spectral_rsl,sigma_psf,wslit)
3233 else if (datatype=='spectrum_euv') then
3234 call get_line_info(spectrum_wl,ion,mass,logte,line_center,spatial_rsl,spectral_rsl,sigma_psf,wslit)
3235 endif
3236
3237 if (mype==0) then
3238 inquire(qunit,opened=fileopen)
3239 if(.not.fileopen)then
3240 ! generate filename
3241 filenr=snapshotini
3242 if (autoconvert) filenr=snapshotnext
3243 if (datatype=='image_euv') then
3244 write(filename,'(a,i4.4,a)') trim(filename_euv),filenr,".vtu"
3245 else if (datatype=='image_sxr') then
3246 write(filename,'(a,i4.4,a)') trim(filename_sxr),filenr,".vtu"
3247 else if (datatype=='image_whitelight') then
3248 write(filename,'(a,i4.4,a)') trim(filename_whitelight),filenr,".vtu"
3249 else if (datatype=='spectrum_euv') then
3250 write(filename,'(a,i4.4,a)') trim(filename_spectrum),filenr,".vtu"
3251 endif
3252 open(qunit,file=filename,status='unknown',form='formatted')
3253 endif
3254 ! generate xml header
3255 write(qunit,'(a)')'<?xml version="1.0"?>'
3256 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
3257 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
3258 write(qunit,'(a)')'<UnstructuredGrid>'
3259 write(qunit,'(a)')'<FieldData>'
3260 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
3261 'NumberOfTuples="1" format="ascii">'
3262 write(qunit,*) real(global_time*time_convert_factor)
3263 write(qunit,'(a)')'</DataArray>'
3264 if (datatype=='image_euv' .or. datatype=='spectrum_euv') then
3265 write(qunit,'(2a)')'<DataArray type="Float32" Name="logT" ',&
3266 'NumberOfTuples="1" format="ascii">'
3267 write(qunit,*) real(logte)
3268 write(qunit,'(a)')'</DataArray>'
3269 endif
3270 write(qunit,'(a)')'</FieldData>'
3271 do ixp=1,npiece
3272 write(qunit,'(a,i7,a,i7,a)') &
3273 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
3274 write(qunit,'(a)')'<CellData>'
3275 do j=1,nwc
3276 if (datatype=='image_euv') then
3277 if (j==1) then
3278 if (wavelength<100) then
3279 write(vname,'(a,i2)') "AIA",wavelength
3280 else if (wavelength<1000) then
3281 write(vname,'(a,i3)') "AIA",wavelength
3282 else
3283 write(vname,'(a,i4)') "IRIS",wavelength
3284 endif
3285 endif
3286 if (j==2 .and. dat_resolution) vname='Doppler_velocity'
3287 else if (datatype=='image_sxr') then
3288 if (emin_sxr<10 .and. emax_sxr<10) then
3289 write(vname,'(a,i1,a,i1,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
3290 else if (emin_sxr<10 .and. emax_sxr>=10) then
3291 write(vname,'(a,i1,a,i2,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
3292 else
3293 write(vname,'(a,i2,a,i2,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
3294 endif
3295 else if (datatype=='image_whitelight') then
3296 write(vname,'(a)')'whitelight'
3297 else if (datatype=='spectrum_euv') then
3298 if (spectrum_wl==1354) then
3299 write(vname,'(a,i4)') "SG",spectrum_wl
3300 else
3301 write(vname,'(a,i3)') "EIS",spectrum_wl
3302 endif
3303 endif
3304 write(qunit,'(a,a,a)')&
3305 '<DataArray type="Float64" Name="',trim(vname),'" format="ascii">'
3306 write(qunit,'(200(1pe14.6))') ((wc(ixp,ixc1,ixc2,j),ixc1=1,nc1),ixc2=1,nc2)
3307 write(qunit,'(a)')'</DataArray>'
3308 enddo
3309 write(qunit,'(a)')'</CellData>'
3310 write(qunit,'(a)')'<Points>'
3311 write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
3312 do ix2=1,np2
3313 do ix1=1,np1
3314 if (datatype=='image_euv' .and. dat_resolution) then
3315 if (los_phi==0 .and. los_theta==90) then
3316 write(qunit,'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
3317 else if (los_phi==90 .and. los_theta==90) then
3318 write(qunit,'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
3319 else
3320 write(qunit,'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
3321 endif
3322 else if (datatype=='image_sxr' .and. dat_resolution) then
3323 if (los_phi==0 .and. los_theta==90) then
3324 write(qunit,'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
3325 else if (los_phi==90 .and. los_theta==90) then
3326 write(qunit,'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
3327 else
3328 write(qunit,'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
3329 endif
3330 else
3331 write(qunit,'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
3332 endif
3333 enddo
3334 enddo
3335 write(qunit,'(a)')'</DataArray>'
3336 write(qunit,'(a)')'</Points>'
3337 ! connetivity part
3338 write(qunit,'(a)')'<Cells>'
3339 write(qunit,'(a)')'<DataArray type="Int32" Name="connectivity" format="ascii">'
3340 do ix2=1,nc2
3341 do ix1=1,nc1
3342 write(qunit,'(4(i7))') ix1-1+(ix2-1)*np1,ix1+(ix2-1)*np1,&
3343 ix1-1+ix2*np1,ix1+ix2*np1
3344 enddo
3345 enddo
3346 write(qunit,'(a)')'</DataArray>'
3347 ! offsets data array
3348 write(qunit,'(a)')'<DataArray type="Int32" Name="offsets" format="ascii">'
3349 do icel=1,nc
3350 write(qunit,'(i7)') icel*(2**2)
3351 enddo
3352 write(qunit,'(a)')'</DataArray>'
3353 ! VTK cell type data array
3354 write(qunit,'(a)')'<DataArray type="Int32" Name="types" format="ascii">'
3355 ! VTK_LINE=3; VTK_PIXEL=8; VTK_VOXEL=11 -> vtk-syntax
3356 vtk_type=8
3357 do icel=1,nc
3358 write(qunit,'(i2)') vtk_type
3359 enddo
3360 write(qunit,'(a)')'</DataArray>'
3361 write(qunit,'(a)')'</Cells>'
3362 write(qunit,'(a)')'</Piece>'
3363 enddo
3364 write(qunit,'(a)')'</UnstructuredGrid>'
3365 write(qunit,'(a)')'</VTKFile>'
3366 close(qunit)
3367 endif
3368 end subroutine write_image_vtucc
3369
3370 subroutine dot_product_loc(vec1,vec2,res)
3371 double precision, intent(in) :: vec1(1:3),vec2(1:3)
3372 double precision, intent(out) :: res
3373
3374 res=vec1(1)*vec2(1)+vec1(2)*vec2(2)+vec1(3)*vec2(3)
3375
3376 end subroutine dot_product_loc
3377
3378 subroutine cross_product_loc(vec_in1,vec_in2,vec_out)
3379 double precision, intent(in) :: vec_in1(1:3),vec_in2(1:3)
3380 double precision, intent(out) :: vec_out(1:3)
3381
3382 vec_out(1)=vec_in1(2)*vec_in2(3)-vec_in1(3)*vec_in2(2)
3383 vec_out(2)=vec_in1(3)*vec_in2(1)-vec_in1(1)*vec_in2(3)
3384 vec_out(3)=vec_in1(1)*vec_in2(2)-vec_in1(2)*vec_in2(1)
3385
3386 end subroutine cross_product_loc
3387
3389 integer :: j
3390 double precision :: LOS_psi
3391 double precision :: vec_car(1:3),vec_z(1:3),vec_temp1(1:3),vec_temp2(1:3)
3392 double precision :: vec_LOS_sph(1:3),vec_xI1_sph(1:3),vec_xI2_sph(1:3)
3393
3394 ! antiparallel to LOS in spherical
3395 vec_los(1)=1.d0
3396 vec_los(2)=dpi*los_theta/180.d0
3397 vec_los(3)=dpi*los_phi/180.d0
3398 ! LOS in cartesian
3399 call spherical_to_cartesian(vec_los,vec_car)
3400 vec_los=-vec_car
3401
3402 ! theta=0 in cartesian
3403 vec_z(:)=zero
3404 vec_z(3)=1.d0
3405
3406 ! x direction for image
3407 if (los_theta==zero) then
3408 vec_temp1(1)=1.d0
3409 vec_temp1(2)=dpi/2.d0
3410 vec_temp1(3)=dpi*los_phi/180.d0
3411 call spherical_to_cartesian(vec_temp1,vec_car)
3412 vec_temp1=-vec_car
3413 call cross_product_loc(vec_temp1,vec_z,vec_xi1)
3414 else
3416 endif
3417
3418 ! y direction for image
3420
3421 ! rotate the image
3422 vec_temp1=vec_xi1/sqrt(vec_xi1(1)**2+vec_xi1(2)**2+vec_xi1(3)**2)
3423 vec_temp2=vec_xi2/sqrt(vec_xi2(1)**2+vec_xi2(2)**2+vec_xi2(3)**2)
3424 los_psi=dpi*image_rotate/180.d0
3425 vec_xi1=vec_temp1*cos(los_psi)-vec_temp2*sin(los_psi)
3426 vec_xi2=vec_temp2*cos(los_psi)+vec_temp1*sin(los_psi)
3427
3428 do j=1,3
3429 if (abs(vec_los(j))<smalldouble) vec_los(j)=zero
3430 if (abs(vec_xi1(j))<smalldouble) vec_xi1(j)=zero
3431 if (abs(vec_xi2(j))<smalldouble) vec_xi2(j)=zero
3432 enddo
3433
3434 call cartesian_to_spherical(vec_los,vec_los_sph)
3435 call cartesian_to_spherical(vec_xi1,vec_xi1_sph)
3436 call cartesian_to_spherical(vec_xi2,vec_xi2_sph)
3437 vec_los_sph(2:3)=vec_los_sph(2:3)*180.d0/dpi
3438 vec_xi1_sph(2:3)=vec_xi1_sph(2:3)*180.d0/dpi
3439 vec_xi2_sph(2:3)=vec_xi2_sph(2:3)*180.d0/dpi
3440
3441 if (mype==0) write(*,'(a,f3.1,f6.1,f6.1,a)') ' LOS vector (spherial): [',vec_los_sph(1),vec_los_sph(2),vec_los_sph(3),']'
3442 if (mype==0) write(*,'(a,f3.1,f6.1,f6.1,a)') ' xI1 vector (spherial): [',vec_xi1_sph(1),vec_xi1_sph(2),vec_xi1_sph(3),']'
3443 if (mype==0) write(*,'(a,f3.1,f6.1,f6.1,a)') ' xI2 vector (spherial): [',vec_xi2_sph(1),vec_xi2_sph(2),vec_xi2_sph(3),']'
3444
3445 end subroutine init_vectors_spherical
3446
3447 subroutine spherical_to_cartesian(vec_sph,vec_car)
3448 ! angles in rad
3449 double precision, intent(in) :: vec_sph(1:3)
3450 double precision, intent(inout) :: vec_car(1:3)
3451
3452 vec_car(1)=vec_sph(1)*dsin(vec_sph(2))*dcos(vec_sph(3))
3453 vec_car(2)=vec_sph(1)*dsin(vec_sph(2))*dsin(vec_sph(3))
3454 vec_car(3)=vec_sph(1)*dcos(vec_sph(2))
3455
3456 end subroutine spherical_to_cartesian
3457
3458 subroutine cartesian_to_spherical(vec_car,vec_sph)
3459 ! angles in rad
3460 double precision, intent(in) :: vec_car(1:3)
3461 double precision, intent(inout) :: vec_sph(1:3)
3462
3463 vec_sph(1)=dsqrt(vec_car(1)**2+vec_car(2)**2+vec_car(3)**2)
3464 vec_sph(2)=dacos(vec_car(3)/vec_sph(1))
3465 vec_sph(3)=atan2(vec_car(2),vec_car(3))
3466
3467 end subroutine cartesian_to_spherical
3468
3470 integer :: j
3471 double precision :: LOS_psi
3472 double precision :: vec_z(1:3),vec_temp1(1:3),vec_temp2(1:3)
3473
3474 ! vectors for image coordinate
3475 vec_los(1)=-cos(dpi*los_phi/180.d0)*sin(dpi*los_theta/180.d0)
3476 vec_los(2)=-sin(dpi*los_phi/180.d0)*sin(dpi*los_theta/180.d0)
3477 vec_los(3)=-cos(dpi*los_theta/180.d0)
3478 do j=1,3
3479 if (abs(vec_los(j))<=smalldouble) vec_los(j)=zero
3480 enddo
3481 vec_z(:)=zero
3482 vec_z(3)=1.d0
3483 if (los_theta==zero) then
3484 vec_xi1(1)=cos(dpi*los_phi/180.d0)
3485 vec_xi1(2)=sin(dpi*los_phi/180.d0)
3486 vec_xi1(3)=zero
3487 else
3489 endif
3491 vec_temp1=vec_xi1/sqrt(vec_xi1(1)**2+vec_xi1(2)**2+vec_xi1(3)**2)
3492 vec_temp2=vec_xi2/sqrt(vec_xi2(1)**2+vec_xi2(2)**2+vec_xi2(3)**2)
3493 los_psi=dpi*image_rotate/180.d0
3494 vec_xi1=vec_temp1*cos(los_psi)-vec_temp2*sin(los_psi)
3495 vec_xi2=vec_temp2*cos(los_psi)+vec_temp1*sin(los_psi)
3496
3497 do j=1,3
3498 if (abs(vec_xi1(j))<smalldouble) vec_xi1(j)=zero
3499 if (abs(vec_xi2(j))<smalldouble) vec_xi2(j)=zero
3500 enddo
3501
3502 if (mype==0) write(*,'(a,f5.2,f6.2,f6.2,a)') ' LOS vector: [',vec_los(1),vec_los(2),vec_los(3),']'
3503 if (mype==0) write(*,'(a,f5.2,f6.2,f6.2,a)') ' xI1 vector: [',vec_xi1(1),vec_xi1(2),vec_xi1(3),']'
3504 if (mype==0) write(*,'(a,f5.2,f6.2,f6.2,a)') ' xI2 vector: [',vec_xi2(1),vec_xi2(2),vec_xi2(3),']'
3505
3506 end subroutine init_vectors_cartesian
3507
3508 subroutine get_cor_image_spherical(x_3D_sph,x_image)
3509 double precision, intent(in) :: x_3D_sph(1:3)
3510 double precision, intent(inout) :: x_image(1:2)
3511 double precision :: res,res_origin
3512 double precision :: x_3D(1:3)
3513
3514 call spherical_to_cartesian(x_3d_sph,x_3d)
3515 call dot_product_loc(x_3d,vec_xi1,res)
3516 x_image(1)=res
3517 call dot_product_loc(x_3d,vec_xi2,res)
3518 x_image(2)=res
3519
3520 end subroutine get_cor_image_spherical
3521
3522 subroutine get_cor_image(x_3D,x_image)
3523 double precision, intent(in) :: x_3D(1:3)
3524 double precision, intent(inout) :: x_image(1:2)
3525 double precision :: res,res_origin
3526
3527 call dot_product_loc(x_3d,vec_xi1,res)
3528 call dot_product_loc(x_origin,vec_xi1,res_origin)
3529 x_image(1)=res-res_origin
3530 call dot_product_loc(x_3d,vec_xi2,res)
3531 call dot_product_loc(x_origin,vec_xi2,res_origin)
3532 x_image(2)=res-res_origin
3533
3534 end subroutine get_cor_image
3535
3536end module mod_thermal_emission
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
double precision, parameter const_rsun
double precision, parameter kb_cgs
Boltzmann constant in cgs.
double precision, parameter half
double precision, parameter one
double precision, parameter dpi
Pi.
double precision, parameter zero
some frequently used numbers
double precision, parameter smalldouble
double precision, parameter mp_cgs
Proton mass in cgs.
double precision, parameter const_c
Module with geometry-related routines (e.g., divergence, curl)
Definition mod_geometry.t:2
integer coordinate
Definition mod_geometry.t:7
integer, parameter spherical
integer, parameter cartesian
Definition mod_geometry.t:8
This module contains definitions of global parameters and variables and some generic functions/subrou...
character(len=std_len) filename_sxr
Base file name for synthetic SXR emission output.
integer spectrum_wl
wave length for spectrum
integer ixghi
Upper index of grid block arrays.
logical activate_unit_arcsec
use arcsec as length unit of images/spectra
character(len=std_len) filename_spectrum
Base file name for synthetic EUV spectrum output.
double precision global_time
The global simulation time.
integer snapshotini
Resume from the snapshot with this index.
character(len=std_len) filename_euv
Base file name for synthetic EUV emission output.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) filename_whitelight
Base file name for synthetic white light.
character(len=std_len) convert_type
Which format to use when converting.
double precision unit_length
Physical scaling factor for length.
double precision location_slit
location of the slit
double precision time_convert_factor
Conversion factor for time unit.
integer icomm
The MPI communicator.
character(len=std_len) whitelight_instrument
white light observation instrument
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
integer snapshotnext
IO: snapshot and collapsed views output numbers/labels.
logical dat_resolution
resolution of the images
double precision r_occultor
the white light emission below it (unit=Rsun) is not visible
integer, dimension(ndim) nstretchedblocks_baselevel
(even) number of (symmetrically) stretched blocks at AMR level 1, per dimension
double precision, dimension(^nd) qstretch_baselevel
stretch factor between cells at AMR level 1, per dimension
double precision unit_velocity
Physical scaling factor for velocity.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision los_theta
direction of the line of sight (LOS)
double precision spectrum_window_max
integer wavelength
wavelength for output
integer, dimension(ndim) stretch_type
What kind of stretching is used per dimension.
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
double precision instrument_resolution_factor
times for enhancing spatial resolution for EUV image/spectra
double precision spectrum_window_min
spectral window
integer refine_max_level
Maximal number of AMR levels.
integer direction_slit
direction of the slit (for dat resolution only)
double precision, dimension(1:3) x_origin
where the is the origin (X=0,Y=0) of image
integer, dimension(:,:), allocatable node
integer, parameter ixglo
Lower index of grid block arrays (always 1)
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition mod_physics.t:4
double precision, dimension(1:3) vec_los
subroutine get_goes_flux_grid(ixil, ixol, w, x, dv, xboxl, xbl, fl, eflux_grid)
subroutine integrate_spectra_cartesian(igrid, wl, dwlg, xs, dxsg, spectra, numwl, numxs, fl)
subroutine get_spectrum_datresol(qunit, datatype, fl)
double precision, dimension(1:101) f_304
double precision, dimension(1:101) f_193
double precision, dimension(1:60) f_264
double precision, dimension(1:60) f_263
subroutine get_euv_image(qunit, fl)
double precision, dimension(1:60) t_eis1
subroutine get_sxr(ixil, ixol, w, x, fl, flux, el, eu)
subroutine get_unit_vector_spherical(x_sph, unitv_r, unitv_theta, unitv_phi)
double precision, dimension(1:101) f_131
subroutine get_line_info(wl, ion, mass, logte, line_center, spatial_px, spectral_px, sigma_psf, width_slit)
double precision, dimension(1:60) f_255
double precision, dimension(1:101) t_aia
subroutine write_image_vtucc(qunit, xc, wc, dxc, npiece, nc1, nc2, nwc, datatype)
subroutine get_cor_image(x_3d, x_image)
subroutine get_thomson_parameters(rl, a, b, c, d)
subroutine integrate_euv_datresol(igrid, nxif1, nxif2, xif1, xif2, dxif1, dxif2, fl, euv, dpl)
double precision, dimension(1:60) t_eis2
double precision, dimension(1:101) f_171
double precision, dimension(1:101) f_94
subroutine integrate_spectra_datresol(igrid, wl, dwl, spectra, numwl, numxs, dir_loc, fl)
double precision, dimension(1:3) vec_xi1
subroutine get_spectrum(qunit, datatype, fl)
subroutine cartesian_to_spherical(vec_car, vec_sph)
subroutine get_cor_image_spherical(x_3d_sph, x_image)
subroutine dot_product_loc(vec1, vec2, res)
subroutine integrate_emission_spherical(igrid, numxi1, numxi2, xi1, xi2, dxi, fl, datatype, em)
subroutine get_image(qunit, datatype, fl)
subroutine integrate_emission_cartesian(igrid, numxi1, numxi2, xi1, xi2, dxi, fl, datatype, em)
subroutine write_image_vticc(qunit, xo1, xo2, dxo1, dxo2, wo, nxo1, nxo2, nwo, nc1, nc2)
subroutine get_sxr_image(qunit, fl)
double precision, dimension(1:60) f_192
subroutine get_whitelight_thomson(rl, rin, ne, a, b, c, d, fluxtb, fluxpb)
subroutine get_image_datresol(qunit, datatype, fl)
subroutine get_goes_sxr_flux(xboxl, fl, eflux)
subroutine integrate_whitelight_spherical(igrid, numxi1, numxi2, numwi, xi1, xi2, dxi, fl, datatype, wlb)
double precision, dimension(1:3) vec_xi2
double precision, dimension(1:41) f_1354
subroutine cross_product_loc(vec_in1, vec_in2, vec_out)
subroutine integrate_sxr_datresol(igrid, nxif1, nxif2, xif1, xif2, dxif1, dxif2, fl, sxr)
double precision, dimension(1:101) f_335
subroutine output_data(qunit, xo1, xo2, dxo1, dxo2, wo, nxo1, nxo2, nwo, datatype)
double precision, dimension(1:41) t_iris
subroutine get_euv_spectrum(qunit, fl)
double precision, dimension(1:101) f_211
subroutine get_whitelight_image(qunit, fl)
subroutine get_euv(wl, ixil, ixol, w, x, fl, flux)
subroutine spherical_to_cartesian(vec_sph, vec_car)