19 double precision ::
t_aia(1:101)
22 double precision ::
f_335(1:101)
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. /
49 data f_94 / 4.25022959
d-37, 4.35880298
d-36, 3.57054296
d-35, 2.18175426
d-34, &
50 8.97592571
d-34, 2.68512961
d-33, 7.49559346
d-33, 2.11603751
d-32, &
51 5.39752853
d-32, 1.02935904
d-31, 1.33822307
d-31, 1.40884290
d-31, &
52 1.54933156
d-31, 2.07543102
d-31, 3.42026227
d-31, 6.31171444
d-31, &
53 1.16559416
d-30, 1.95360497
d-30, 2.77818735
d-30, 3.43552578
d-30, &
54 4.04061803
d-30, 4.75470982
d-30, 5.65553769
d-30, 6.70595782
d-30, &
55 7.80680354
d-30, 8.93247715
d-30, 1.02618156
d-29, 1.25979030
d-29, &
56 1.88526483
d-29, 3.62448572
d-29, 7.50553279
d-29, 1.42337571
d-28, &
57 2.37912813
d-28, 3.55232305
d-28, 4.84985757
d-28, 6.20662827
d-28, &
58 7.66193687
d-28, 9.30403645
d-28, 1.10519802
d-27, 1.25786927
d-27, &
59 1.34362634
d-27, 1.33185242
d-27, 1.22302081
d-27, 1.05677973
d-27, &
60 9.23064720
d-28, 8.78570994
d-28, 8.02397416
d-28, 5.87681142
d-28, &
61 3.82272695
d-28, 3.11492649
d-28, 3.85736090
d-28, 5.98893519
d-28, &
62 9.57553548
d-28, 1.46650267
d-27, 2.10365847
d-27, 2.79406671
d-27, &
63 3.39420087
d-27, 3.71077520
d-27, 3.57296767
d-27, 2.95114380
d-27, &
64 2.02913103
d-27, 1.13361825
d-27, 5.13405629
d-28, 2.01305089
d-28, &
65 8.15781482
d-29, 4.28366817
d-29, 3.08701543
d-29, 2.68693906
d-29, &
66 2.51764203
d-29, 2.41773103
d-29, 2.33996083
d-29, 2.26997246
d-29, &
67 2.20316143
d-29, 2.13810001
d-29, 2.07424438
d-29, 2.01149189
d-29, &
68 1.94980213
d-29, 1.88917920
d-29, 1.82963583
d-29, 1.77116920
d-29, &
69 1.71374392
d-29, 1.65740593
d-29, 1.60214447
d-29, 1.54803205
d-29, &
70 1.49510777
d-29, 1.44346818
d-29, 1.39322305
d-29, 1.34441897
d-29, &
71 1.29713709
d-29, 1.25132618
d-29, 1.20686068
d-29, 1.14226584
d-29, &
72 1.09866413
d-29, 1.05635524
d-29, 1.01532444
d-29, 9.75577134
d-30, &
73 9.37102736
d-30, 8.99873335
d-30, 8.63860172
d-30, 8.29051944
d-30, &
76 data f_131 / 3.18403601
d-37, 3.22254703
d-36, 2.61657920
d-35, &
77 1.59575286
d-34, 6.65779556
d-34, 2.07015132
d-33, &
78 6.05768615
d-33, 1.76074833
d-32, 4.52633001
d-32, &
79 8.57121883
d-32, 1.09184271
d-31, 1.10207963
d-31, &
80 1.11371658
d-31, 1.29105226
d-31, 1.80385897
d-31, &
81 3.27295431
d-31, 8.92002136
d-31, 3.15214579
d-30, &
82 9.73440787
d-30, 2.22709702
d-29, 4.01788984
d-29, &
83 6.27471832
d-29, 8.91764995
d-29, 1.18725647
d-28, &
84 1.52888040
d-28, 2.05082946
d-28, 3.47651873
d-28, &
85 8.80482184
d-28, 2.66533063
d-27, 7.05805149
d-27, &
86 1.46072515
d-26, 2.45282476
d-26, 3.55303726
d-26, &
87 4.59075911
d-26, 5.36503515
d-26, 5.68444094
d-26, &
88 5.47222296
d-26, 4.81119761
d-26, 3.85959059
d-26, &
89 2.80383406
d-26, 1.83977650
d-26, 1.11182849
d-26, &
90 6.50748885
d-27, 3.96843481
d-27, 2.61876319
d-27, &
91 1.85525324
d-27, 1.39717024
d-27, 1.11504283
d-27, &
92 9.38169611
d-28, 8.24801234
d-28, 7.43331919
d-28, &
93 6.74537063
d-28, 6.14495760
d-28, 5.70805277
d-28, &
94 5.61219786
d-28, 6.31981777
d-28, 9.19747307
d-28, &
95 1.76795732
d-27, 3.77985446
d-27, 7.43166191
d-27, &
96 1.19785603
d-26, 1.48234676
d-26, 1.36673114
d-26, &
97 9.61047146
d-27, 5.61209353
d-27, 3.04779780
d-27, &
98 1.69378976
d-27, 1.02113491
d-27, 6.82223774
d-28, &
99 5.02099099
d-28, 3.99377760
d-28, 3.36279037
d-28, &
100 2.94767378
d-28, 2.65740865
d-28, 2.44396277
d-28, &
101 2.28003967
d-28, 2.14941419
d-28, 2.04178995
d-28, &
102 1.95031045
d-28, 1.87011994
d-28, 1.79777869
d-28, &
103 1.73093957
d-28, 1.66795789
d-28, 1.60785455
d-28, &
104 1.55002399
d-28, 1.49418229
d-28, 1.44022426
d-28, &
105 1.38807103
d-28, 1.33772767
d-28, 1.28908404
d-28, &
106 1.24196208
d-28, 1.17437501
d-28, 1.12854330
d-28, &
107 1.08410498
d-28, 1.04112003
d-28, 9.99529904
d-29, &
108 9.59358806
d-29, 9.20512291
d-29, 8.83009123
d-29, &
109 8.46817043
d-29, 8.11921928
d-29 /
111 data f_171 / 2.98015581
d-42, 1.24696230
d-40, 3.37614652
d-39, 5.64103034
d-38, &
112 5.20550266
d-37, 2.77785939
d-36, 1.16283616
d-35, 6.50007689
d-35, &
113 9.96177399
d-34, 1.89586076
d-32, 2.10982799
d-31, 1.36946479
d-30, &
114 6.27396553
d-30, 2.29955134
d-29, 7.13430211
d-29, 1.91024282
d-28, &
115 4.35358848
d-28, 7.94807808
d-28, 1.07431875
d-27, 1.08399488
d-27, &
116 9.16212938
d-28, 7.34715770
d-28, 6.59246382
d-28, 9.13541375
d-28, &
117 2.05939035
d-27, 5.08206555
d-27, 1.10148083
d-26, 2.01884662
d-26, &
118 3.13578384
d-26, 4.14367719
d-26, 5.36067711
d-26, 8.74170213
d-26, &
119 1.64161233
d-25, 2.94587860
d-25, 4.76298332
d-25, 6.91765639
d-25, &
120 9.08825111
d-25, 1.08496183
d-24, 1.17440114
d-24, 1.13943939
d-24, &
121 9.71696981
d-25, 7.09593688
d-25, 4.31376399
d-25, 2.12708486
d-25, &
122 8.47429567
d-26, 3.17608104
d-26, 1.95898842
d-26, 1.98064242
d-26, &
123 1.67706555
d-26, 8.99126003
d-27, 3.29773878
d-27, 1.28896127
d-27, &
124 8.51169698
d-28, 7.53520167
d-28, 6.18268143
d-28, 4.30034650
d-28, &
125 2.78152409
d-28, 1.95437088
d-28, 1.65896278
d-28, 1.68740181
d-28, &
126 1.76054383
d-28, 1.63978419
d-28, 1.32880591
d-28, 1.00833205
d-28, &
127 7.82252806
d-29, 6.36181741
d-29, 5.34633869
d-29, 4.58013864
d-29, &
128 3.97833422
d-29, 3.49414760
d-29, 3.09790940
d-29, 2.76786227
d-29, &
129 2.48806269
d-29, 2.24823367
d-29, 2.04016653
d-29, 1.85977413
d-29, &
130 1.70367499
d-29, 1.56966125
d-29, 1.45570643
d-29, 1.35964565
d-29, &
131 1.27879263
d-29, 1.21016980
d-29, 1.15132499
d-29, 1.09959628
d-29, &
132 1.05307482
d-29, 1.01040261
d-29, 9.70657096
d-30, 9.33214234
d-30, &
133 8.97689427
d-30, 8.63761192
d-30, 8.31149879
d-30, 7.85162401
d-30, &
134 7.53828281
d-30, 7.23559452
d-30, 6.94341530
d-30, 6.66137038
d-30, &
135 6.38929156
d-30, 6.12669083
d-30, 5.87346434
d-30, 5.62943622
d-30, &
138 data f_193 / 6.40066486
d-32, 4.92737300
d-31, 2.95342934
d-30, 1.28061594
d-29, &
139 3.47747667
d-29, 5.88554792
d-29, 7.72171179
d-29, 9.75609282
d-29, &
140 1.34318963
d-28, 1.96252638
d-28, 2.70163878
d-28, 3.63192965
d-28, &
141 5.28087341
d-28, 8.37821446
d-28, 1.39089159
d-27, 2.31749718
d-27, &
142 3.77510689
d-27, 5.85198594
d-27, 8.26021568
d-27, 1.04870405
d-26, &
143 1.25209374
d-26, 1.47406787
d-26, 1.77174067
d-26, 2.24098537
d-26, &
144 3.05926105
d-26, 4.50018853
d-26, 6.84720216
d-26, 1.00595861
d-25, &
145 1.30759222
d-25, 1.36481773
d-25, 1.15943558
d-25, 1.01467304
d-25, &
146 1.04092532
d-25, 1.15071251
d-25, 1.27416033
d-25, 1.38463476
d-25, &
147 1.47882726
d-25, 1.57041238
d-25, 1.69786224
d-25, 1.94970397
d-25, &
148 2.50332918
d-25, 3.58321431
d-25, 5.18061550
d-25, 6.60405549
d-25, &
149 6.64085365
d-25, 4.83825816
d-25, 2.40545020
d-25, 8.59534098
d-26, &
150 2.90920638
d-26, 1.33204845
d-26, 9.03933926
d-27, 7.78910836
d-27, &
151 7.29342321
d-27, 7.40267022
d-27, 8.05279981
d-27, 8.13829291
d-27, &
152 6.92634262
d-27, 5.12521880
d-27, 3.59527615
d-27, 2.69617560
d-27, &
153 2.84432713
d-27, 5.06697306
d-27, 1.01281903
d-26, 1.63526978
d-26, &
154 2.06759342
d-26, 2.19482312
d-26, 2.10050611
d-26, 1.89837248
d-26, &
155 1.66347131
d-26, 1.43071097
d-26, 1.21518419
d-26, 1.02078343
d-26, &
156 8.46936184
d-27, 6.93015742
d-27, 5.56973237
d-27, 4.38951754
d-27, &
157 3.38456457
d-27, 2.55309556
d-27, 1.88904224
d-27, 1.38057546
d-27, &
158 1.00718330
d-27, 7.43581116
d-28, 5.63562931
d-28, 4.43359435
d-28, &
159 3.63923535
d-28, 3.11248143
d-28, 2.75586846
d-28, 2.50672237
d-28, &
160 2.32419348
d-28, 2.18325682
d-28, 2.06834486
d-28, 1.93497044
d-28, &
161 1.84540751
d-28, 1.76356504
d-28, 1.68741425
d-28, 1.61566157
d-28, &
162 1.54754523
d-28, 1.48249410
d-28, 1.42020176
d-28, 1.36045230
d-28, &
165 data f_211 / 4.74439912
d-42, 1.95251522
d-40, 5.19700194
d-39, 8.53120166
d-38, &
166 7.72745727
d-37, 4.04158559
d-36, 1.64853511
d-35, 8.56295439
d-35, &
167 1.17529722
d-33, 2.16867729
d-32, 2.40472264
d-31, 1.56418133
d-30, &
168 7.20032889
d-30, 2.65838271
d-29, 8.33196904
d-29, 2.26128236
d-28, &
169 5.24295811
d-28, 9.77791121
d-28, 1.35913489
d-27, 1.43957785
d-27, &
170 1.37591544
d-27, 1.49029886
d-27, 2.06183401
d-27, 3.31440622
d-27, &
171 5.42497318
d-27, 8.41100374
d-27, 1.17941366
d-26, 1.49269794
d-26, &
172 1.71506074
d-26, 1.71266353
d-26, 1.51434781
d-26, 1.36766622
d-26, &
173 1.33483562
d-26, 1.36834518
d-26, 1.45829002
d-26, 1.62575306
d-26, &
174 1.88773347
d-26, 2.22026986
d-26, 2.54930499
d-26, 2.80758138
d-26, &
175 3.06176409
d-26, 3.62799792
d-26, 5.13226109
d-26, 8.46260744
d-26, &
176 1.38486586
d-25, 1.86192535
d-25, 1.78007934
d-25, 1.16548409
d-25, &
177 5.89293257
d-26, 2.69952884
d-26, 1.24891081
d-26, 6.41273176
d-27, &
178 4.08282914
d-27, 3.26463328
d-27, 2.76230280
d-27, 2.08986882
d-27, &
179 1.37658470
d-27, 8.48489381
d-28, 5.19304217
d-28, 3.19312514
d-28, &
180 2.02968197
d-28, 1.50171666
d-28, 1.39164218
d-28, 1.42448821
d-28, &
181 1.41714519
d-28, 1.33341059
d-28, 1.20759270
d-28, 1.07259692
d-28, &
182 9.44895400
d-29, 8.29030041
d-29, 7.25440631
d-29, 6.33479483
d-29, &
183 5.51563757
d-29, 4.79002469
d-29, 4.14990482
d-29, 3.59384972
d-29, &
184 3.12010860
d-29, 2.72624742
d-29, 2.40734791
d-29, 2.15543565
d-29, &
185 1.95921688
d-29, 1.80682882
d-29, 1.68695662
d-29, 1.59020936
d-29, &
186 1.50940886
d-29, 1.43956179
d-29, 1.37731622
d-29, 1.32049043
d-29, &
187 1.26771875
d-29, 1.21803879
d-29, 1.17074716
d-29, 1.10507836
d-29, &
188 1.06022834
d-29, 1.01703080
d-29, 9.75436986
d-30, 9.35349257
d-30, &
189 8.96744546
d-30, 8.59527489
d-30, 8.23678940
d-30, 7.89144480
d-30, &
192 data f_304 / 3.62695850
d-32, 2.79969087
d-31, 1.68340584
d-30, 7.32681440
d-30, &
193 1.99967770
d-29, 3.41296785
d-29, 4.55409104
d-29, 5.94994635
d-29, &
194 8.59864963
d-29, 1.39787633
d-28, 3.17701965
d-28, 1.14474920
d-27, &
195 4.44845958
d-27, 1.54785841
d-26, 4.70265345
d-26, 1.24524365
d-25, &
196 2.81535352
d-25, 5.10093666
d-25, 6.83545307
d-25, 6.82110329
d-25, &
197 5.66886188
d-25, 4.36205513
d-25, 3.29265688
d-25, 2.49802368
d-25, &
198 1.92527113
d-25, 1.51058572
d-25, 1.20596047
d-25, 9.76884267
d-26, &
199 7.89979266
d-26, 6.18224289
d-26, 4.67298332
d-26, 3.57934505
d-26, &
200 2.84535785
d-26, 2.32853022
d-26, 1.95228514
d-26, 1.67880071
d-26, &
201 1.47608785
d-26, 1.32199691
d-26, 1.20070960
d-26, 1.09378177
d-26, &
202 1.00031730
d-26, 9.62434001
d-27, 1.05063954
d-26, 1.27267143
d-26, &
203 1.45923057
d-26, 1.36746707
d-26, 1.03466970
d-26, 6.97647829
d-27, &
204 4.63141039
d-27, 3.19031994
d-27, 2.33373613
d-27, 1.81589079
d-27, &
205 1.48446917
d-27, 1.26611478
d-27, 1.12617468
d-27, 1.03625148
d-27, &
206 9.61400595
d-28, 8.79016231
d-28, 7.82612130
d-28, 6.73762960
d-28, &
207 5.59717956
d-28, 4.53010243
d-28, 3.65712196
d-28, 3.00958686
d-28, &
208 2.54011502
d-28, 2.18102277
d-28, 1.88736437
d-28, 1.63817539
d-28, &
209 1.42283147
d-28, 1.23631916
d-28, 1.07526003
d-28, 9.36797928
d-29, &
210 8.18565660
d-29, 7.18152734
d-29, 6.32523238
d-29, 5.59513985
d-29, &
211 4.96614048
d-29, 4.42518826
d-29, 3.95487628
d-29, 3.54690294
d-29, &
212 3.18953930
d-29, 2.87720933
d-29, 2.60186750
d-29, 2.36011522
d-29, &
213 2.14717806
d-29, 1.95905217
d-29, 1.79287981
d-29, 1.64562262
d-29, &
214 1.51489425
d-29, 1.39876064
d-29, 1.29496850
d-29, 1.18665438
d-29, &
215 1.10240474
d-29, 1.02643099
d-29, 9.57780996
d-30, 8.95465151
d-30, &
216 8.38950190
d-30, 7.87283711
d-30, 7.40136507
d-30, 6.96804279
d-30, &
219 data f_335 / 2.46882661
d-32, 1.89476632
d-31, 1.13216502
d-30, 4.89532008
d-30, &
220 1.32745970
d-29, 2.25390335
d-29, 3.00511672
d-29, 3.96035934
d-29, &
221 5.77977656
d-29, 8.58600736
d-29, 1.14083000
d-28, 1.48644411
d-28, &
222 2.15788823
d-28, 3.51628877
d-28, 6.12200698
d-28, 1.08184987
d-27, &
223 1.85590697
d-27, 2.91679107
d-27, 3.94405396
d-27, 4.63610680
d-27, &
224 5.13824456
d-27, 5.66602209
d-27, 6.30009232
d-27, 7.03422868
d-27, &
225 7.77973918
d-27, 8.32371831
d-27, 8.56724316
d-27, 8.62601374
d-27, &
226 8.13308844
d-27, 6.53188216
d-27, 4.55197029
d-27, 3.57590087
d-27, &
227 3.59571707
d-27, 4.03502770
d-27, 4.54366411
d-27, 4.96914990
d-27, &
228 5.24601170
d-27, 5.39979250
d-27, 5.43023669
d-27, 5.26235042
d-27, &
229 4.91585495
d-27, 4.52628362
d-27, 4.13385020
d-27, 3.67538967
d-27, &
230 3.39939742
d-27, 3.81284533
d-27, 5.02332701
d-27, 6.19438602
d-27, &
231 6.49613071
d-27, 6.04010475
d-27, 5.24664275
d-27, 4.37225997
d-27, &
232 3.52957182
d-27, 2.76212276
d-27, 2.08473158
d-27, 1.50850518
d-27, &
233 1.04602472
d-27, 7.13091243
d-28, 5.34289645
d-28, 5.21079581
d-28, &
234 6.22246365
d-28, 6.99555864
d-28, 6.29665489
d-28, 4.45077026
d-28, &
235 2.67046793
d-28, 1.52774686
d-28, 9.18061770
d-29, 6.09116074
d-29, &
236 4.48562572
d-29, 3.59463696
d-29, 3.05820218
d-29, 2.70766652
d-29, &
237 2.46144034
d-29, 2.27758450
d-29, 2.13331183
d-29, 2.01537836
d-29, &
238 1.91566180
d-29, 1.82893912
d-29, 1.75167748
d-29, 1.68136168
d-29, &
239 1.61615595
d-29, 1.55481846
d-29, 1.49643236
d-29, 1.44046656
d-29, &
240 1.38657085
d-29, 1.33459068
d-29, 1.28447380
d-29, 1.23615682
d-29, &
241 1.18963296
d-29, 1.14478976
d-29, 1.10146637
d-29, 1.04039479
d-29, &
242 9.98611410
d-30, 9.58205147
d-30, 9.19202009
d-30, 8.81551313
d-30, &
243 8.45252127
d-30, 8.10224764
d-30, 7.76469090
d-30, 7.43954323
d-30, &
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, &
259 data f_1354 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
260 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
261 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
262 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
263 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
264 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 1.09503647
d-39, &
265 5.47214550
d-36, 2.42433983
d-33, 2.75295034
d-31, 1.21929718
d-29, &
266 2.48392125
d-28, 2.33268145
d-27, 8.68623633
d-27, 1.00166284
d-26, &
267 3.63126633
d-27, 7.45174807
d-28, 1.38224064
d-28, 2.69270994
d-29, &
268 5.53314977
d-30, 1.15313092
d-30, 2.34195788
d-31, 4.48242942
d-32, &
274 data t_eis1 / 1.99526231
d+05, 2.23872114
d+05, 2.51188643
d+05, 2.81838293
d+05, &
275 3.16227766
d+05, 3.54813389
d+05, 3.98107171
d+05, 4.46683592
d+05, &
276 5.01187234
d+05, 5.62341325
d+05, 6.30957344
d+05, 7.07945784
d+05, &
277 7.94328235
d+05, 8.91250938
d+05, 1.00000000
d+06, 1.12201845
d+06, &
278 1.25892541
d+06, 1.41253754
d+06, 1.58489319
d+06, 1.77827941
d+06, &
279 1.99526231
d+06, 2.23872114
d+06, 2.51188643
d+06, 2.81838293
d+06, &
280 3.16227766
d+06, 3.54813389
d+06, 3.98107171
d+06, 4.46683592
d+06, &
281 5.01187234
d+06, 5.62341325
d+06, 6.30957344
d+06, 7.07945784
d+06, &
282 7.94328235
d+06, 8.91250938
d+06, 1.00000000
d+07, 1.12201845
d+07, &
283 1.25892541
d+07, 1.41253754
d+07, 1.58489319
d+07, 1.77827941
d+07, &
284 1.99526231
d+07, 2.23872114
d+07, 2.51188643
d+07, 2.81838293
d+07, &
285 3.16227766
d+07, 3.54813389
d+07, 3.98107171
d+07, 4.46683592
d+07, &
286 5.01187234
d+07, 5.62341325
d+07, 6.30957344
d+07, 7.07945784
d+07, &
287 7.94328235
d+07, 8.91250938
d+07, 1.00000000
d+08, 1.12201845
d+08, &
288 1.25892541
d+08, 1.41253754
d+08, 1.58489319
d+08, 1.77827941
d+08 /
290 data t_eis2 / 1.99526231
d+06, 2.23872114
d+06, 2.51188643
d+06, 2.81838293
d+06, &
291 3.16227766
d+06, 3.54813389
d+06, 3.98107171
d+06, 4.46683592
d+06, &
292 5.01187234
d+06, 5.62341325
d+06, 6.30957344
d+06, 7.07945784
d+06, &
293 7.94328235
d+06, 8.91250938
d+06, 1.00000000
d+07, 1.12201845
d+07, &
294 1.25892541
d+07, 1.41253754
d+07, 1.58489319
d+07, 1.77827941
d+07, &
295 1.99526231
d+07, 2.23872114
d+07, 2.51188643
d+07, 2.81838293
d+07, &
296 3.16227766
d+07, 3.54813389
d+07, 3.98107171
d+07, 4.46683592
d+07, &
297 5.01187234
d+07, 5.62341325
d+07, 6.30957344
d+07, 7.07945784
d+07, &
298 7.94328235
d+07, 8.91250938
d+07, 1.00000000
d+08, 1.12201845
d+08, &
299 1.25892541
d+08, 1.41253754
d+08, 1.58489319
d+08, 1.77827941
d+08, &
300 1.99526231
d+08, 2.23872114
d+08, 2.51188643
d+08, 2.81838293
d+08, &
301 3.16227766
d+08, 3.54813389
d+08, 3.98107171
d+08, 4.46683592
d+08, &
302 5.01187234
d+08, 5.62341325
d+08, 6.30957344
d+08, 7.07945784
d+08, &
303 7.94328235
d+08, 8.91250938
d+08, 1.00000000
d+09, 1.12201845
d+09, &
304 1.25892541
d+09, 1.41253754
d+09, 1.58489319
d+09, 1.77827941
d+09 /
306 data f_263 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
307 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
308 0.00000000
d+00, 4.46454917
d-45, 3.26774829
d-42, 1.25292566
d-39, &
309 2.66922338
d-37, 3.28497742
d-35, 2.38677554
d-33, 1.03937729
d-31, &
310 2.75075687
d-30, 4.47961733
d-29, 4.46729177
d-28, 2.64862689
d-27, &
311 8.90863800
d-27, 1.72437548
d-26, 2.22217752
d-26, 2.27999477
d-26, &
312 2.08264363
d-26, 1.78226687
d-26, 1.45821699
d-26, 1.14675379
d-26, &
313 8.63082492
d-27, 6.15925429
d-27, 4.11252514
d-27, 2.51530564
d-27, &
314 1.37090986
d-27, 6.42443134
d-28, 2.48392636
d-28, 7.59187874
d-29, &
315 1.77852938
d-29, 3.23945221
d-30, 4.90533903
d-31, 6.75458158
d-32, &
316 9.06878868
d-33, 1.23927474
d-33, 1.75769395
d-34, 2.60710914
d-35, &
317 4.04318030
d-36, 6.53500581
d-37, 1.09365022
d-37, 1.88383322
d-38, &
318 3.31425233
d-39, 5.90964084
d-40, 1.06147549
d-40, 1.90706170
d-41, &
319 3.41331584
d-42, 6.07310718
d-43, 1.07364738
d-43, 1.89085498
d-44, &
320 3.32598922
d-45, 5.87125640
d-46, 0.00000000
d+00, 0.00000000
d+00 /
322 data f_264 / 0.00000000
d+00, 2.81670057
d-46, 1.28007268
d-43, 2.54586603
d-41, &
323 2.67887256
d-39, 1.68413285
d-37, 6.85702304
d-36, 1.91797284
d-34, &
324 3.84675839
d-33, 5.69939170
d-32, 6.36224608
d-31, 5.39176489
d-30, &
325 3.45478458
d-29, 1.64848693
d-28, 5.71476364
d-28, 1.39909997
d-27, &
326 2.37743056
d-27, 2.86712530
d-27, 2.65206348
d-27, 2.07175767
d-27, &
327 1.47866767
d-27, 1.01087374
d-27, 6.79605811
d-28, 4.54746770
d-28, &
328 3.04351751
d-28, 2.03639149
d-28, 1.35940991
d-28, 9.01451939
d-29, &
329 5.91289972
d-29, 3.81821178
d-29, 2.41434696
d-29, 1.48871220
d-29, &
330 8.93362094
d-30, 5.21097445
d-30, 2.95964719
d-30, 1.64278748
d-30, &
331 8.95571660
d-31, 4.82096011
d-31, 2.57390991
d-31, 1.36821781
d-31, &
332 7.27136350
d-32, 3.87019426
d-32, 2.06883430
d-32, 1.11228884
d-32, &
333 6.01883313
d-33, 3.27790676
d-33, 1.79805012
d-33, 9.93085346
d-34, &
334 5.52139556
d-34, 3.08881387
d-34, 1.73890315
d-34, 9.84434964
d-35, &
335 5.60603378
d-35, 3.20626492
d-35, 1.84111068
d-35, 0.00000000
d+00, &
336 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00 /
338 data f_192 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 4.35772105
d-44, &
339 1.26162319
d-41, 1.97471205
d-39, 1.83409019
d-37, 1.08206288
d-35, &
340 4.27914363
d-34, 1.17943846
d-32, 2.32565755
d-31, 3.33087991
d-30, &
341 3.47013260
d-29, 2.60375866
d-28, 1.37737127
d-27, 5.01053913
d-27, &
342 1.23479810
d-26, 2.11310542
d-26, 2.71831513
d-26, 2.89851163
d-26, &
343 2.77312376
d-26, 2.50025229
d-26, 2.18323661
d-26, 1.86980322
d-26, &
344 1.58035034
d-26, 1.31985651
d-26, 1.08733133
d-26, 8.81804906
d-27, &
345 7.00417973
d-27, 5.43356567
d-27, 4.09857884
d-27, 2.99651764
d-27, &
346 2.11902962
d-27, 1.45014127
d-27, 9.62291023
d-28, 6.21548647
d-28, &
347 3.92807578
d-28, 2.44230375
d-28, 1.50167782
d-28, 9.17611405
d-29, &
348 5.58707641
d-29, 3.40570915
d-29, 2.08030862
d-29, 1.27588676
d-29, &
349 7.86535588
d-30, 4.87646151
d-30, 3.03888897
d-30, 1.90578649
d-30, &
350 1.20195947
d-30, 7.61955060
d-31, 4.85602199
d-31, 3.11049969
d-31, &
351 2.00087065
d-31, 1.29223740
d-31, 8.37422008
d-32, 0.00000000
d+00, &
352 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00 /
354 data f_255 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 1.76014287
d-44, &
355 5.07057938
d-42, 7.90473970
d-40, 7.31852999
d-38, 4.30709255
d-36, &
356 1.70009061
d-34, 4.67925160
d-33, 9.21703546
d-32, 1.31918676
d-30, &
357 1.37393161
d-29, 1.03102379
d-28, 5.45694018
d-28, 1.98699648
d-27, &
358 4.90346776
d-27, 8.40524725
d-27, 1.08321456
d-26, 1.15714525
d-26, &
359 1.10905152
d-26, 1.00155023
d-26, 8.75799161
d-27, 7.50935839
d-27, &
360 6.35253533
d-27, 5.30919268
d-27, 4.37669455
d-27, 3.55185164
d-27, &
361 2.82347055
d-27, 2.19257595
d-27, 1.65589541
d-27, 1.21224987
d-27, &
362 8.58395132
d-28, 5.88163935
d-28, 3.90721447
d-28, 2.52593407
d-28, &
363 1.59739995
d-28, 9.93802874
d-29, 6.11343388
d-29, 3.73711135
d-29, &
364 2.27618743
d-29, 1.38793199
d-29, 8.48060787
d-30, 5.20305940
d-30, &
365 3.20867365
d-30, 1.99011277
d-30, 1.24064551
d-30, 7.78310544
d-31, &
366 4.91013681
d-31, 3.11338381
d-31, 1.98451675
d-31, 1.27135460
d-31, &
367 8.17917486
d-32, 5.28280497
d-32, 3.42357159
d-32, 0.00000000
d+00, &
368 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00 /
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)
384 integer,
intent(in) :: ixI^L, ixO^L
385 double precision,
intent(in) :: w(ixI^S, nw)
386 double precision,
intent(in) :: x(ixI^S, 1:ndim)
387 double precision,
intent(out):: val1(ixI^S), val2(ixI^S)
393 procedure(
get_subr1),
pointer,
nopass :: get_rho => null()
394 procedure(
get_subr1),
pointer,
nopass :: get_pthermal => null()
395 procedure(
get_subr1),
pointer,
nopass :: get_var_rfactor => null()
406 double precision,
allocatable :: source(:^d&)
407 double precision,
allocatable :: opacity(:^d&)
408 double precision,
allocatable :: sourcev(:^d&)
409 double precision,
allocatable :: xface1(:),xface2(:),xface3(:)
410 double precision,
allocatable :: rface(:),thetaface(:),phiface(:)
411 double precision,
allocatable :: rface2(:),theta_cos(:),phi_sin(:),phi_cos(:)
412 double precision :: box_min(1:3)=0.d0
413 double precision :: box_max(1:3)=0.d0
418 logical :: has_pixels=.false.
430 character(len=*),
intent(in) :: datatype
433 call mpistop(
"bad radiation_transfer")
438 call mpistop(
"dat_resolution_mode must be nominal or minimum")
454 case(
'cart',
'cart_dda')
456 case(
'spherical',
'sph_intersection')
458 case(
'sph_dda',
'spherical_dda')
468 call mpistop(
"bad emission_model")
472 call mpistop(
"tau and absorption-fraction output need thick transfer")
476 call mpistop(
"radsyn_pixel_batch must be positive")
479 call mpistop(
"radsyn_segment_batch_factor must be non-negative")
482 call mpistop(
"radsyn_segment_memory_mb must be positive")
485 call mpistop(
"radsyn_segment_comm_factor must be positive")
490 call mpistop(
"instrument_postprocess currently needs dat-resolution EUV images")
493 call mpistop(
"instrument_postprocess is not yet supported for spherical rays")
496 call mpistop(
"instrument_postprocess currently supports only EUV AIA or radio_ff images")
499 call mpistop(
"radio_ff instrument_postprocess needs radio_beam_fwhm > 0 arcsec")
507 if (datatype /=
'image_euv' .and. datatype /=
'spectrum_euv')
then
508 call mpistop(
"emission_model=euv_aia is only valid for EUV synthesis")
511 if (datatype /=
'image_whitelight')
then
512 call mpistop(
"emission_model=white_light is only valid for white-light synthesis")
515 if (datatype /=
'image_euv')
then
516 call mpistop(
"emission_model=radio_ff currently reuses EUV-image convert types")
519 call mpistop(
"emission_model=radio_ff needs radio_frequency > 0")
521 case(
'pseudo_current')
522 if (datatype /=
'image_euv')
then
523 call mpistop(
"emission_model=pseudo_current is only valid for EUV-image convert types")
526 call mpistop(
"emission_model=pseudo_current currently supports only thin transfer")
531 if (datatype /=
'image_euv' .or. .not.
slab)
then
532 call mpistop(
"ray_method=cart needs Cartesian EUV slab images")
537 call mpistop(
"ray_method=spherical currently needs 3D spherical grids")
540 call mpistop(
"ray_method=spherical currently needs 3D spherical grids")
545 call mpistop(
"bad ray_method=spherical mode")
548 call mpistop(
"ray_method=spherical currently supports only EUV AIA emission")
550 if (xprobmin2<=1.
d-10 .or. xprobmax2>=dpi-1.
d-10)
then
551 call mpistop(
"ray_method=spherical does not support polar-axis crossing domains")
553 if (xprobmax3<=xprobmin3 .or. xprobmax3-xprobmin3>=2.d0*dpi-1.
d-10)
then
554 call mpistop(
"ray_method=spherical does not support phi-wrapping domains")
559 if (datatype /=
'image_euv')
then
560 call mpistop(
"thick transfer is only defined for EUV images")
565 if (.not.
slab)
call mpistop(
"cartesian thick EUV currently needs slab output")
567 call mpistop(
"thick EUV currently needs Cartesian dat_resolution output")
573 call mpistop(
"thick EUV currently needs x/y/z-aligned LOS")
580 double precision,
intent(in) :: emissivity,opacity,path_length
581 double precision,
intent(inout) :: intensity,tau
583 double precision :: dtau
585 if (path_length<=zero)
return
587 dtau=max(zero,opacity)*path_length
593 trim(emission_model)/=
'radio_ff' .and. &
598 logical,
intent(in) :: has_doppler,has_thick
601 if (has_doppler) num_outputs=num_outputs+1
602 if (has_thick .and. output_tau) num_outputs=num_outputs+1
603 if (has_thick .and. output_absorption_fraction) num_outputs=num_outputs+1
607 integer,
intent(in) :: nI1,nI2
608 double precision,
intent(in) :: EUV(nI1,nI2),unitv
609 double precision,
intent(inout) :: Dpl(nI1,nI2)
615 if (euv(ix1,ix2)/=zero)
then
616 dpl(ix1,ix2)=(dpl(ix1,ix2)/euv(ix1,ix2))*unitv
620 if (abs(dpl(ix1,ix2))<smalldouble) dpl(ix1,ix2)=zero
626 integer,
intent(in) :: nI1,nI2
627 double precision,
intent(in) :: EUV(nI1,nI2),EUVthin(nI1,nI2),smallflux
628 double precision,
intent(out) :: Absorption(nI1,nI2)
629 logical,
intent(in),
optional :: cap_to_one
632 logical :: cap_absorption
635 cap_absorption=.false.
636 if (
present(cap_to_one)) cap_absorption=cap_to_one
639 if (euvthin(ix1,ix2)>smallflux)
then
640 absorption(ix1,ix2)=max(zero,(euvthin(ix1,ix2)-euv(ix1,ix2))/euvthin(ix1,ix2))
641 if (cap_absorption) absorption(ix1,ix2)=min(one,absorption(ix1,ix2))
647 subroutine pack_euv_image_outputs(nI1,nI2,EUV,wI,smallflux,has_doppler,has_thick,Dpl,Tau,EUVthin,&
649 integer,
intent(in) :: nI1,nI2
650 double precision,
intent(in) :: EUV(nI1,nI2),smallflux
651 double precision,
intent(inout) :: wI(:,:,:)
652 logical,
intent(in) :: has_doppler,has_thick
653 double precision,
intent(in),
optional :: Dpl(nI1,nI2),Tau(nI1,nI2),EUVthin(nI1,nI2)
654 logical,
intent(in),
optional :: cap_absorption
657 double precision,
allocatable :: Absorption(:,:)
662 if (has_doppler)
then
663 if (.not.
present(dpl))
call mpistop(
"Doppler output requested without Doppler image")
667 if (has_thick .and. output_tau)
then
668 if (.not.
present(tau))
call mpistop(
"tau output requested without tau image")
672 if (has_thick .and. output_absorption_fraction)
then
673 if (.not.
present(euvthin))
call mpistop(
"absorption output requested without thin image")
674 allocate(absorption(ni1,ni2))
677 wi(:,:,iw)=absorption(:,:)
678 deallocate(absorption)
683 integer,
intent(out) :: pixel_batch_target,segment_batch_target,segment_comm_target
685 pixel_batch_target=max(1,radsyn_pixel_batch)
686 if (radsyn_segment_batch_factor>0)
then
687 segment_batch_target=max(128,radsyn_segment_batch_factor*pixel_batch_target)
689 segment_batch_target=max(128,int(min(dble(huge(segment_batch_target)),&
690 max(128.d0,radsyn_segment_memory_mb*1048576.d0/256.d0))))
692 segment_comm_target=max(128,radsyn_segment_comm_factor*pixel_batch_target)
696 double precision,
intent(in) :: tau
706 double precision,
intent(in) :: argument
708 if (argument<-700.d0)
then
710 else if (argument>700.d0)
then
718 double precision,
intent(in) :: exponent
720 if (exponent>300.d0)
then
722 else if (exponent<-300.d0)
then
730 double precision,
intent(in) :: temperature
731 integer,
intent(in) :: n_table
732 double precision,
intent(in) :: t_table(n_table),f_table(n_table)
733 logical,
intent(in) :: log_temperature,log_response
735 integer :: ilo,ihi,imid
736 double precision :: temp_lookup,response_lookup,flo,fhi
739 if (temperature<=zero)
return
740 if (log_temperature)
then
741 temp_lookup=log10(temperature)
743 temp_lookup=temperature
745 if (temp_lookup<t_table(1) .or. temp_lookup>t_table(n_table))
return
746 if (temp_lookup==t_table(n_table))
then
747 if (log_response)
then
748 response_lookup=log10(max(f_table(n_table),1.d-99))
750 response_lookup=f_table(n_table)
757 if (temp_lookup>=t_table(imid))
then
763 if (log_response)
then
764 flo=log10(max(f_table(ilo),1.d-99))
765 fhi=log10(max(f_table(ilo+1),1.d-99))
770 response_lookup=flo*(temp_lookup-t_table(ilo+1))/(t_table(ilo)-t_table(ilo+1))+&
771 fhi*(temp_lookup-t_table(ilo))/(t_table(ilo+1)-t_table(ilo))
774 if (log_response)
then
783 integer,
intent(in) :: ixI^L, ixO^L, n_table
784 double precision,
intent(in) :: Te(ixI^S),t_table(n_table),f_table(n_table)
785 double precision,
intent(inout) :: flux(ixI^S)
786 logical,
intent(in) :: log_temperature,log_response
789 double precision :: GT
791 {
do ix^db=ixomin^db,ixomax^db\}
793 flux(ix^d)=flux(ix^d)*gt
794 if (flux(ix^d)<zero) flux(ix^d)=zero
804 double precision,
intent(in) :: Te,Ne
805 double precision,
intent(out) :: x_HII,x_HeII,x_HeIII
807 double precision :: Pe,log_H21,log_He21,log_He32,log_He321
808 double precision :: logScaleHe,w_H21,term0,term1,term2,denHe
809 double precision,
parameter :: Xe_H21=13.6d0
810 double precision,
parameter :: Xe_He21=24.587d0
811 double precision,
parameter :: Xe_He32=54.416d0
821 log_h21=2.5d0*log10(te)-5040.d0*xe_h21/te-log10(pe)-0.48d0
822 log_he21=log10(4.d0)+2.5d0*log10(te)-5040.d0*xe_he21/te-log10(pe)-0.48d0
823 log_he32=2.5d0*log10(te)-5040.d0*xe_he32/te-log10(pe)-0.48d0
826 x_hii=w_h21/(1.d0+w_h21)
830 log_he321=log_he21+log_he32
831 logscalehe=max(
zero,log_he21,log_he321)
835 denhe=term0+term1+term2
848 double precision,
intent(in) :: nH,Te,rHe,Ne_guess
849 double precision,
intent(out) :: x_HII,x_HeII,x_HeIII
851 integer,
parameter :: max_iter=32
853 double precision :: Ne,Ne_lo,Ne_hi,Ne_new,residual,derivative
854 double precision :: e_He,de_HII_dNe,de_He_dNe
859 if (nh<=zero .or. te<=zero)
return
861 ne_lo=max(1.d-30*nh,1.d-100)
862 ne_hi=(1.d0+2.d0*rhe)*nh
863 ne=min(max(ne_guess,ne_lo),ne_hi)
867 e_he=x_heii+2.d0*x_heiii
868 residual=ne/nh-x_hii-rhe*e_he
869 if (abs(residual)<1.d-10)
exit
871 if (residual>zero)
then
879 de_hii_dne=-x_hii*(1.d0-x_hii)/ne
880 de_he_dne=(x_heii*(e_he-1.d0) &
881 +2.d0*x_heiii*(e_he-2.d0))/ne
882 derivative=1.d0/nh-de_hii_dne-rhe*de_he_dne
883 ne_new=ne-residual/derivative
884 if (.not.(ne_new>ne_lo .and. ne_new<ne_hi))
then
885 ne_new=0.5d0*(ne_lo+ne_hi)
899 integer,
intent(in) :: wl
900 integer,
intent(in) :: ixI^L, ixO^L
901 double precision,
intent(in) :: x(ixI^S,1:ndim)
902 double precision,
intent(in) :: w(ixI^S,1:nw)
904 double precision,
intent(out) :: kappa(ixI^S)
907 double precision :: pth(ixI^S),rho(ixI^S),Rfactor(ixI^S),Te(ixI^S)
908 double precision :: Ne(ixI^S),nH(ixI^S)
909 double precision :: wave_ratio,s_H1,s_He1,s_He2
910 double precision :: x_HII,x_HeII,x_HeIII,iz_H,iz_He,Rdummy
911 double precision :: N_H1,N_He1,N_He2
912 double precision,
parameter :: rHe_opacity=0.1d0
913 double precision,
parameter :: sigma_H1=5.16d-20, sigma_he1=9.25d-19, sigma_he2=7.17d-19
915 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
916 call fl%get_rho(w,x,ixi^l,ixo^l,rho)
917 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,rfactor)
919 {
do ix^db=ixomin^db,ixomax^db\}
920 if (rho(ix^d)>zero .and. rfactor(ix^d)>zero)
then
921 te(ix^d)=pth(ix^d)/(rho(ix^d)*rfactor(ix^d))*unit_temperature
928 call eos%get_ne_nH(ixi^l, ixo^l, w, x, ne, nh)
930 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
931 nh(ixo^s)=nh(ixo^s)*unit_numberdensity/1.d6
933 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
934 nh(ixo^s)=nh(ixo^s)*unit_numberdensity
937 wave_ratio=dble(wl)/171.d0
941 if (wl<=912) s_h1=wave_ratio**3*sigma_h1
942 if (wl<=504) s_he1=wave_ratio**2*sigma_he1
943 if (wl<=228) s_he2=wave_ratio**2.75d0*sigma_he2
946 {
do ix^db=ixomin^db,ixomax^db\}
947 if (te(ix^d)>zero .and. nh(ix^d)>zero)
then
948 select case (trim(eos%eos_type))
956 x_heii=iz_he*(1.d0-iz_he)
963 if (ne(ix^d)>zero)
then
965 x_hii,x_heii,x_heiii)
966 if (trim(eos%method)==
'analytic')
then
969 x_hii=min(one,max(zero,ne(ix^d)/nh(ix^d)))
971 x_hii=min(one,max(zero,ne(ix^d)/nh(ix^d) &
972 -eos%He_abundance*(x_heii+2.d0*x_heiii)))
976 rhe_opacity,ne(ix^d),x_hii,x_heii,x_heiii)
983 rhe_opacity,ne(ix^d),x_hii,x_heii,x_heiii)
986 n_h1=nh(ix^d)*(1.d0-x_hii)
987 n_he1=rhe_opacity*nh(ix^d)*(1.d0-x_heii-x_heiii)
988 n_he2=rhe_opacity*nh(ix^d)*x_heii
989 kappa(ix^d)=max(zero,n_h1*s_h1+n_he1*s_he1+n_he2*s_he2)
995 integer,
intent(in) :: igrid
996 integer,
intent(in) :: ixI^L, ixO^L
997 double precision,
intent(in) :: w(ixI^S,1:nw)
998 double precision,
intent(out) :: source(ixI^S)
1000 integer :: ix^D,idir,idirmin,idirmin0
1001 double precision :: current(ixI^S,7-2*ndir:3)
1003 if (.not.
allocated(iw_mag))
then
1004 call mpistop(
"emission_model=pseudo_current needs magnetic-field variables")
1009 call curlvector(w(ixi^s,iw_mag(1:ndir)),ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
1011 current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+ps(igrid)%J0(ixo^s,idirmin0:3)
1015 {
do ix^db=ixomin^db,ixomax^db\}
1017 source(ix^d)=source(ix^d)+current(ix^d,idir)**2
1025 integer,
intent(in) :: ixI^L, ixO^L
1026 double precision,
intent(in) :: x(ixI^S,1:ndim)
1027 double precision,
intent(in) :: w(ixI^S,1:nw)
1029 double precision,
intent(out) :: source(ixI^S),kappa(ixI^S)
1032 double precision :: pth(ixI^S),Te(ixI^S),Ne(ixI^S)
1033 double precision :: nH_dummy(ixI^S),gff
1035 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
1036 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
1037 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
1038 te(ixo^s)=pth(ixo^s)/(ne(ixo^s)*te(ixo^s))*unit_temperature
1039 call eos%get_ne_nH(ixi^l,ixo^l,w, x,ne,nh_dummy)
1041 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
1043 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
1048 {
do ix^db=ixomin^db,ixomax^db\}
1049 if (te(ix^d)>zero .and. ne(ix^d)>zero)
then
1050 if (te(ix^d)<2.d5)
then
1051 gff=18.2d0+1.5d0*log(te(ix^d))-log(radio_frequency)
1053 gff=24.5d0+log(te(ix^d))-log(radio_frequency)
1056 kappa(ix^d)=9.78d-3*ne(ix^d)**2*gff/(radio_frequency**2*te(ix^d)**1.5d0)
1057 source(ix^d)=te(ix^d)*kappa(ix^d)
1062 subroutine get_line_info(wl,ion,mass,logTe,line_center,spatial_px,spectral_px,sigma_PSF,width_slit)
1074 integer,
intent(in) :: wl
1075 integer,
intent(out) :: mass
1076 character(len=30),
intent(out) :: ion
1077 double precision,
intent(out) :: logTe,line_center,spatial_px,spectral_px
1078 double precision,
intent(out) :: sigma_PSF,width_slit
1148 line_center=1354.1d0
1150 spectral_px=12.98
d-3
1157 line_center=262.976d0
1166 line_center=263.765d0
1175 line_center=192.028d0
1184 line_center=255.113d0
1190 call mpistop(
"No information about this line")
1204 integer,
intent(in) :: wl
1205 integer,
intent(in) :: ixI^L, ixO^L
1206 double precision,
intent(in) :: x(ixI^S,1:ndim)
1207 double precision,
intent(in) :: w(ixI^S,1:nw)
1209 double precision,
intent(out) :: flux(ixI^S)
1212 double precision :: pth(ixI^S),rho(ixI^S),Rfactor(ixI^S),Te(ixI^S)
1213 double precision :: Ne(ixI^S),nH(ixI^S)
1215 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
1216 call fl%get_rho(w,x,ixi^l,ixo^l,rho)
1217 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,rfactor)
1223 call eos%get_ne_nH(ixi^l, ixo^l, w, x, ne, nh)
1231 flux(ixo^s)=ne(ixo^s)*nh(ixo^s)
1259 call mpistop(
"Unknown wavelength")
1271 integer,
intent(in) :: ixI^L,ixO^L
1272 integer,
intent(in) :: El,Eu
1273 double precision,
intent(in) :: x(ixI^S,1:ndim)
1274 double precision,
intent(in) :: w(ixI^S,nw)
1276 double precision,
intent(out) :: flux(ixI^S)
1278 integer :: ix^D,ixO^D
1280 double precision :: I0,kb,keV,dE,Ei
1281 double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
1282 double precision :: Ne(ixI^S),gff(ixI^S),fi(ixI^S)
1283 double precision :: EM(ixI^S)
1284 double precision :: nH_dummy(ixI^S)
1290 nume=floor((eu-el)/de)
1291 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
1292 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
1293 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
1296 call eos%get_ne_nH(ixi^l, ixo^l, w, x, ne, nh_dummy)
1299 em(ixo^s)=(ne(ixo^s))**2*1.d6
1302 em(ixo^s)=(ne(ixo^s))**2
1304 kbt(ixo^s)=kb*te(ixo^s)/kev
1309 {
do ix^db=ixomin^db,ixomax^db\}
1310 if (kbt(ix^d)>0.01*ei)
then
1311 if(kbt(ix^d)<ei) gff(ix^d)=(kbt(ix^d)/ei)**0.4
1312 fi(ix^d)=(em(ix^d)*gff(ix^d))*
exp_clamped(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
1317 flux(ixo^s)=flux(ixo^s)+fi(ixo^s)*de
1319 flux(ixo^s)=flux(ixo^s)*i0
1326 double precision,
intent(in) :: xbox^L
1328 double precision,
intent(out) :: eflux
1330 double precision :: dxb^D,xb^L
1331 integer :: iigrid,igrid,j
1332 integer :: ixO^L,ixI^L,ix^D
1333 double precision :: eflux_grid,eflux_pe
1335 ^d&iximin^d=
ixglo^d;
1336 ^d&iximax^d=
ixghi^d;
1337 ^d&ixomin^d=ixmlo^d;
1338 ^d&ixomax^d=ixmhi^d;
1340 do iigrid=1,igridstail; igrid=igrids(iigrid);
1344 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)
1345 eflux_pe=eflux_pe+eflux_grid
1347 call mpi_allreduce(eflux_pe,eflux,1,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
1354 integer,
intent(in) :: ixI^L,ixO^L
1355 double precision,
intent(in) :: x(ixI^S,1:ndim),dV(ixI^S)
1356 double precision,
intent(in) :: w(ixI^S,nw)
1357 double precision,
intent(in) :: xbox^L,xb^L
1359 double precision,
intent(out) :: eflux_grid
1361 integer :: ix^D,ixO^D,ixb^L
1362 integer :: iE,numE,j,inbox
1363 double precision :: I0,kb,keV,dE,Ei,El,Eu,A_cgs
1364 double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
1365 double precision :: Ne(ixI^S),EM(ixI^S)
1366 double precision :: nH_dummy(ixI^S)
1367 double precision :: gff,fi,erg_SI
1371 {
if (xbmin^d<xboxmax^d .and. xbmax^d>xboxmin^d) inbox=inbox+1\}
1373 if (inbox==ndim)
then
1375 ^d&ixbmin^d=ixomin^d;
1376 ^d&ixbmax^d=ixomax^d;
1377 {
if (xbmax^d>xboxmax^d) ixbmax^d=ixomax^d-ceiling((xbmax^d-xboxmax^d)/
dxlevel(^d))\}
1378 {
if (xbmin^d<xboxmin^d) ixbmin^d=ceiling((xboxmin^d-xbmin^d)/
dxlevel(^d))+ixomin^d\}
1385 el=const_h*const_c/(8.d0*a_cgs)/kev
1386 eu=const_h*const_c/(1.d0*a_cgs)/kev
1388 nume=floor((eu-el)/de)
1389 call fl%get_pthermal(w,x,ixi^l,ixb^l,pth)
1390 call fl%get_rho(w,x,ixi^l,ixb^l,ne)
1391 call fl%get_var_Rfactor(w,x,ixi^l,ixb^l,te)
1394 call eos%get_ne_nH(ixi^l, ixb^l, w, x, ne, nh_dummy)
1397 em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*(
unit_length*1.d2)**3
1400 em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*
unit_length**3
1402 kbt(ixb^s)=kb*te(ixb^s)/kev
1407 {
do ix^db=ixbmin^db,ixbmax^db\}
1408 if (kbt(ix^d)>1.d-2*ei)
then
1409 if(kbt(ix^d)<ei)
then
1410 gff=(kbt(ix^d)/ei)**0.4
1414 fi=(em(ix^d)*gff)*
exp_clamped(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
1415 eflux_grid=eflux_grid+fi*de*ei
1419 eflux_grid=eflux_grid*kev*erg_si
1428 integer,
intent(in) :: qunit
1430 character(20) :: datatype
1433 character (30) :: ion
1434 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1435 double precision :: xslit,arcsec
1437 datatype=
'spectrum_euv'
1442 if (
mype==0) print *,
'###################################################'
1445 if (
mype==0) print *,
'Systhesizing EUV spectrum (observed by IRIS).'
1446 case (263,264,192,255)
1447 if (
mype==0) print *,
'Systhesizing EUV spectrum (observed by Hinode/EIS).'
1449 call mpistop(
'Wrong wavelength!')
1453 call mpistop(
'Wrong spectrum window!')
1456 if (
mype==0)
write(*,
'(a,f8.3,a)')
' Wavelength: ',linecent,
' Angstrom'
1457 if (
mype==0) print *,
'Unit of EUV flux: DN s^-1 pixel^-1'
1461 write(*,
'(a,f5.3,a,f5.1,a)')
' Supposed pixel: ',wlrsl,
' Angstrom x ',spacersl*725.0,
' km'
1462 print *,
'Unit of wavelength: Angstrom (0.1 nm) '
1464 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
1466 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
1468 write(*,
'(a,f8.1,a)')
' Supposed width of slit: ',wslit*725.0,
' km'
1473 print *,
'Unit of wavelength: Angstrom (0.1 nm) '
1475 write(*,
'(a,f5.3,a,f5.1,a)')
' Pixel: ',wlrsl,
' Angstrom x ',spacersl*725.0,
' km'
1476 print *,
'Unit of length: arcsec (~725 km)'
1477 write(*,
'(a,f8.1,a)')
' Location of slit: xI1 = ',
location_slit,
' arcsec'
1478 write(*,
'(a,f8.1,a)')
' Width of slit: ',wslit,
' arcsec'
1481 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
1483 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
1485 write(*,
'(a,f8.1,a)')
' Location of slit: xI1 = ',
location_slit,
' Unit_length'
1486 write(*,
'(a,f8.1,a)')
' Width of slit: ',wslit*725.d0,
' km'
1489 if (
mype==0) print *,
'Direction of the slit: parallel to xI2 vector'
1493 call mpistop(
"EUV spectrum synthesis: support for sperical coordinates is to be added!")
1497 if (
mype==0) print *,
'###################################################'
1503 integer,
intent(in) :: qunit
1504 character(20),
intent(in) :: datatype
1507 integer :: numWL,numXS,iwL,ixS,numWI,numS
1508 double precision :: dwLg,xSmin,xSmax,wLmin,wLmax
1509 double precision,
allocatable :: wL(:),xS(:),dwL(:),dxS(:)
1510 double precision,
allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
1511 integer :: strtype,nstrb,nbb,nuni,nstr,bnx
1512 double precision :: qs,dxfirst,dxmid,lenstr
1514 integer :: iigrid,igrid,j,dir_loc
1515 double precision :: xbmin(1:ndim),xbmax(1:ndim)
1518 numwl=4*int((spectrum_window_max-spectrum_window_min)/(4.d0*dwlg))
1519 wlmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwlg*numwl/2
1520 wlmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwlg*numwl/2
1521 allocate(wl(numwl),dwl(numwl))
1524 wl(iwl)=wlmin+iwl*dwlg-half*dwlg
1527 select case(direction_slit)
1529 numxs=domain_nx1*2**(refine_max_level-1)
1534 strtype=stretch_type(1)
1535 nstrb=nstretchedblocks_baselevel(1)
1536 qs=qstretch_baselevel(1)
1537 if (mype==0) print *,
'Direction of the slit: x'
1539 numxs=domain_nx2*2**(refine_max_level-1)
1544 strtype=stretch_type(2)
1545 nstrb=nstretchedblocks_baselevel(2)
1546 qs=qstretch_baselevel(2)
1547 if (mype==0) print *,
'Direction of the slit: y'
1549 numxs=domain_nx3*2**(refine_max_level-1)
1554 strtype=stretch_type(3)
1555 nstrb=nstretchedblocks_baselevel(3)
1556 qs=qstretch_baselevel(3)
1557 if (mype==0) print *,
'Direction of the slit: z'
1559 call mpistop(
'Wrong direction_slit')
1562 allocate(xs(numxs),dxs(numxs),spectra(numwl,numxs),spectra_rc(numwl,numxs))
1564 allocate(wi(numwl,numxs,numwi))
1566 select case(strtype)
1568 dxs(:)=(xsmax-xsmin)/numxs
1570 xs(ixs)=xsmin+dxs(ixs)*(ixs-half)
1573 qs=qs**(one/2**(refine_max_level-1))
1574 dxfirst=(xsmax-xsmin)*(one-qs)/(one-qs**numxs)
1577 dxs(ixs)=dxfirst*qs**(ixs-1)
1578 xs(ixs)=dxs(1)/(one-qs)*(one-qs**(ixs-1))+half*dxs(ixs)
1584 lenstr=(xsmax-xsmin)/(2.d0+nuni*(one-qs)/(one-qs**nstr))
1585 dxfirst=(xsmax-xsmin)/(dble(nuni)+2.d0/(one-qs)*(one-qs**nstr))
1588 nstr=nstr*2**(refine_max_level-1)
1589 nuni=nuni*2**(refine_max_level-1)
1590 qs=qs**(one/2**(refine_max_level-1))
1591 dxfirst=lenstr*(one-qs)/(one-qs**nstr)
1592 dxmid=dxmid/2**(refine_max_level-1)
1594 if(nuni .gt. 0)
then
1595 do ixs=nstr+1,nstr+nuni
1597 xs(ixs)=lenstr+(dble(ixs)-0.5d0-nstr)*dxs(ixs)+xsmin
1602 dxs(ixs)=dxfirst*qs**(nstr-ixs)
1603 xs(ixs)=xsmin+lenstr-dxs(ixs)*half-dxfirst*(one-qs**(nstr-ixs))/(one-qs)
1606 do ixs=nstr+nuni+1,numxs
1607 dxs(ixs)=dxfirst*qs**(ixs-nstr-nuni-1)
1608 xs(ixs)=xsmax-lenstr+dxs(ixs)*half+dxfirst*(one-qs**(ixs-nstr-nuni-1))/(one-qs)
1611 call mpistop(
"unknown stretch type")
1614 if (los_phi==0 .and. los_theta==90 .and. direction_slit==2)
then
1617 else if (los_phi==0 .and. los_theta==90 .and. direction_slit==3)
then
1620 else if (los_phi==90 .and. los_theta==90 .and. direction_slit==1)
then
1623 else if (los_phi==90 .and. los_theta==90 .and. direction_slit==3)
then
1626 else if (los_theta==0 .and. direction_slit==1)
then
1629 else if (los_theta==0 .and. direction_slit==2)
then
1633 call mpistop(
'Wrong combination of LOS and slit direction!')
1636 if (dir_loc==1)
then
1637 if (location_slit>xprobmax1 .or. location_slit<xprobmin1)
then
1638 call mpistop(
'Wrong value for location_slit!')
1640 if(mype==0)
write(*,
'(a,f8.1,a)')
' Location of slit: x = ',location_slit,
' Unit_length'
1641 else if (dir_loc==2)
then
1642 if (location_slit>xprobmax2 .or. location_slit<xprobmin2)
then
1643 call mpistop(
'Wrong value for location_slit!')
1645 if(mype==0)
write(*,
'(a,f8.1,a)')
' Location of slit: y = ',location_slit,
' Unit_length'
1647 if (location_slit>xprobmax3 .or. location_slit<xprobmin3)
then
1648 call mpistop(
'Wrong value for location_slit!')
1650 if(mype==0)
write(*,
'(a,f8.1,a)')
' Location of slit: z = ',location_slit,
' Unit_length'
1655 do iigrid=1,igridstail; igrid=igrids(iigrid);
1656 ^d&xbmin(^d)=rnode(rpxmin^d_,igrid);
1657 ^d&xbmax(^d)=rnode(rpxmax^d_,igrid);
1658 if (location_slit>=xbmin(dir_loc) .and. location_slit<xbmax(dir_loc))
then
1664 call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
1665 mpi_sum,icomm,ierrmpi)
1668 if (spectra_rc(iwl,ixs)>smalldouble)
then
1669 wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
1676 call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
1678 deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
1685 integer,
intent(in) :: igrid,numWL,numXS,dir_loc
1687 double precision,
intent(in) :: wL(numWL),dwL(numWL)
1688 double precision,
intent(inout) :: spectra(numWL,numXS)
1690 integer :: direction_LOS
1691 integer :: ixO^L,ixI^L,ix^D,ixOnew
1692 double precision,
allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
1693 double precision :: wlc,wlwd
1696 double precision :: logTe,lineCent
1697 character (30) :: ion
1698 double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
1700 integer :: levelg,rft,ixSmin,ixSmax,iwL
1701 double precision :: flux_pix,dL
1703 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1705 if (los_phi==0 .and. los_theta==90)
then
1707 else if (los_phi==90 .and. los_theta==90)
then
1713 ^d&ixomin^d=ixmlo^d\
1714 ^d&ixomax^d=ixmhi^d\
1715 ^d&iximin^d=ixglo^d\
1716 ^d&iximax^d=ixghi^d\
1717 allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
1720 if (dir_loc==1)
then
1721 do ix1=ixomin1,ixomax1
1722 if (location_slit>=(ps(igrid)%x(ix^d,1)-
half*ps(igrid)%dx(ix^d,1)) .and. &
1723 location_slit<(ps(igrid)%x(ix^d,1)+
half*ps(igrid)%dx(ix^d,1)))
then
1729 else if (dir_loc==2)
then
1730 do ix2=ixomin2,ixomax2
1731 if (location_slit>=(ps(igrid)%x(ix^d,2)-
half*ps(igrid)%dx(ix^d,2)) .and. &
1732 location_slit<(ps(igrid)%x(ix^d,2)+
half*ps(igrid)%dx(ix^d,2)))
then
1739 do ix3=ixomin3,ixomax3
1740 if (location_slit>=(ps(igrid)%x(ix^d,3)-
half*ps(igrid)%dx(ix^d,3)) .and. &
1741 location_slit<(ps(igrid)%x(ix^d,3)+
half*ps(igrid)%dx(ix^d,3)))
then
1749 call get_euv(spectrum_wl,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
1750 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
1751 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
1752 v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
1753 call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
1754 call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
1755 te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
1758 levelg=ps(igrid)%level
1759 rft=2**(refine_max_level-levelg)
1761 {
do ix^d=ixomin^d,ixomax^d\}
1764 wlc=linecent*(1.d0+v(ix^d)*unit_velocity*1.d2/
const_c)
1766 wlc=linecent*(1.d0+v(ix^d)*unit_velocity/
const_c)
1768 wlwd=sqrt(
kb_cgs*te(ix^d)*unit_temperature/(mass*
mp_cgs))
1771 select case(direction_slit)
1773 ixsmin=(block_nx1*(node(pig1_,igrid)-1)+(ix1-ixomin1))*rft+1
1774 ixsmax=(block_nx1*(node(pig1_,igrid)-1)+(ix1-ixomin1+1))*rft
1776 ixsmin=(block_nx2*(node(pig2_,igrid)-1)+(ix2-ixomin2))*rft+1
1777 ixsmax=(block_nx2*(node(pig2_,igrid)-1)+(ix2-ixomin2+1))*rft
1779 ixsmin=(block_nx3*(node(pig3_,igrid)-1)+(ix3-ixomin3))*rft+1
1780 ixsmax=(block_nx3*(node(pig3_,igrid)-1)+(ix3-ixomin3+1))*rft
1783 select case(direction_los)
1785 dl=ps(igrid)%dx(ix^d,1)*unit_length
1787 dl=ps(igrid)%dx(ix^d,2)*unit_length
1789 dl=ps(igrid)%dx(ix^d,3)*unit_length
1791 if (si_unit) dl=dl*1.d2
1794 flux_pix=flux(ix^d)*wlrsl*dl*
exp_clamped(-(wl(iwl)-wlc)**2/(2*wlwd**2))/(sqrt(2*
dpi)*wlwd)
1796 flux_pix=flux_pix*wslit/spacersl
1797 spectra(iwl,ixsmin:ixsmax)=spectra(iwl,ixsmin:ixsmax)+flux_pix
1803 deallocate(flux,v,pth,te,rho)
1809 integer,
intent(in) :: qunit
1810 character(20),
intent(in) :: datatype
1813 integer :: numWL,numXS,iwL,ixS,numWI,ix^D
1814 double precision :: dwLg,dxSg,xSmin,xSmax,xScent,wLmin,wLmax
1815 double precision,
allocatable :: wL(:),xS(:),dwL(:),dxS(:)
1816 double precision,
allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
1817 double precision :: vec_cor(1:3),xI_cor(1:2)
1818 double precision :: res,r_loc,r_max
1821 character (30) :: ion
1822 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1823 double precision :: unitv,arcsec,RHESSI_rsl,pixel
1824 integer :: iigrid,igrid,i,j,numS
1825 double precision :: xLmin,xLmax,xslit
1836 xsmin=-abs(xprobmax1)
1837 xsmax=abs(xprobmax1)
1840 if (ix1==1) vec_cor(1)=xprobmin1
1841 if (ix1==2) vec_cor(1)=xprobmax1
1843 if (ix2==1) vec_cor(2)=xprobmin2
1844 if (ix2==2) vec_cor(2)=xprobmax2
1846 if (ix3==1) vec_cor(3)=xprobmin3
1847 if (ix3==2) vec_cor(3)=xprobmax3
1849 r_loc=(vec_cor(1)-x_origin(1))**2
1850 r_loc=r_loc+(vec_cor(2)-x_origin(2))**2
1851 r_loc=r_loc+(vec_cor(3)-x_origin(3))**2
1853 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1856 r_max=max(r_max,r_loc)
1860 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1864 xsmin=min(xsmin,xi_cor(2))
1865 xsmax=max(xsmax,xi_cor(2))
1876 xscent=(xsmin+xsmax)/2.d0
1880 arcsec=7.25d5/unit_length
1882 arcsec=7.25d7/unit_length
1884 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1885 dxsg=spacersl*arcsec
1886 numxs=ceiling((xsmax-xscent)/dxsg)
1887 xsmin=xscent-numxs*dxsg
1888 xsmax=xscent+numxs*dxsg
1891 numwl=2*int((spectrum_window_max-spectrum_window_min)/(2.d0*dwlg))
1892 wlmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwlg*numwl/2
1893 wlmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwlg*numwl/2
1894 allocate(wl(numwl),dwl(numwl),xs(numxs),dxs(numxs))
1896 allocate(wi(numwl,numxs,numwi),spectra(numwl,numxs),spectra_rc(numwl,numxs))
1898 wl(iwl)=wlmin+iwl*dwlg-half*dwlg
1902 xs(ixs)=xsmin+dxsg*(ixs-half)
1908 do iigrid=1,igridstail; igrid=igrids(iigrid);
1910 if (ix1==1) vec_cor(1)=rnode(rpxmin1_,igrid)
1911 if (ix1==2) vec_cor(1)=rnode(rpxmax1_,igrid)
1913 if (ix2==1) vec_cor(2)=rnode(rpxmin2_,igrid)
1914 if (ix2==2) vec_cor(2)=rnode(rpxmax2_,igrid)
1916 if (ix3==1) vec_cor(3)=rnode(rpxmin3_,igrid)
1917 if (ix3==2) vec_cor(3)=rnode(rpxmax3_,igrid)
1919 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1923 xlmin=min(xlmin,xi_cor(1))
1924 xlmax=max(xlmax,xi_cor(1))
1930 if (activate_unit_arcsec)
then
1931 xslit=location_slit*arcsec
1935 if (xslit>=xlmin-wslit*arcsec .and. xslit<=xlmax+wslit*arcsec)
then
1941 call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
1942 mpi_sum,icomm,ierrmpi)
1945 if (spectra_rc(iwl,ixs)>smalldouble)
then
1946 wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
1953 if (activate_unit_arcsec)
then
1958 call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
1960 deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
1966 integer,
intent(in) :: igrid,numWL,numXS
1967 double precision,
intent(in) :: wL(numWL),xS(numXS)
1968 double precision,
intent(in) :: dwLg,dxSg
1969 double precision,
intent(inout) :: spectra(numWL,numXS)
1972 integer :: ixO^L,ixI^L,ix^D,ixOnew,j
1973 double precision,
allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
1974 double precision :: wlc,wlwd,res,dst_slit,xslit,arcsec
1975 double precision :: vloc(1:3),xloc(1:3),dxloc(1:3),xIloc(1:2),dxIloc(1:2)
1976 integer :: nSubC^D,iSubC^D,iwL,ixS,ixSmin,ixSmax,iwLmin,iwLmax,nwL
1977 double precision :: slit_width,dxSubC^D,xerf^L,fluxSubC
1978 double precision :: xSubC(1:3),xCent(1:2)
1981 double precision :: logTe,lineCent
1982 character (30) :: ion
1983 double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
1984 double precision :: sigma_wl,sigma_xs,factor
1987 arcsec=7.25d5/unit_length
1989 arcsec=7.25d7/unit_length
1991 if (activate_unit_arcsec)
then
1992 xslit=location_slit*arcsec
1997 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1999 ^d&ixomin^d=ixmlo^d\
2000 ^d&ixomax^d=ixmhi^d\
2001 ^d&iximin^d=ixglo^d\
2002 ^d&iximax^d=ixghi^d\
2003 allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
2005 call get_euv(spectrum_wl,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
2006 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
2007 call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
2008 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
2009 call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
2010 te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
2011 {
do ix^d=ixomin^d,ixomax^d\}
2013 vloc(j)=ps(igrid)%w(ix^d,iw_mom(j))/rho(ix^d)
2021 slit_width=wslit*arcsec
2022 sigma_wl=sigma_psf*dwlg
2023 sigma_xs=sigma_psf*dxsg
2024 {
do ix^d=ixomin^d,ixomax^d\}
2025 if (flux(ix^d)>smalldouble)
then
2026 xloc(1:3)=ps(igrid)%x(ix^d,1:3)
2027 dxloc(1:3)=ps(igrid)%dx(ix^d,1:3)
2031 if (xiloc(1)>=xslit-half*(slit_width+dxiloc(1)) .and. &
2032 xiloc(1)<=xslit+half*(slit_width+dxiloc(1)))
then
2034 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi1(^d))/(slit_width/16.d0)));
2035 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi2(^d))/(dxsg/4.d0)));
2036 ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
2039 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length*1.d2/dxsg/dxsg
2040 wlc=linecent*(1.d0+v(ix^d)*unit_velocity*1.d2/const_c)
2042 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length/dxsg/dxsg
2043 wlc=linecent*(1.d0+v(ix^d)*unit_velocity/const_c)
2045 wlwd=sqrt(kb_cgs*te(ix^d)*unit_temperature/(mass*mp_cgs))
2046 wlwd=wlwd*linecent/const_c
2048 {
do isubc^d=1,nsubc^d\}
2049 ^d&xsubc(^d)=xloc(^d)-half*dxloc(^d)+(isubc^d-half)*dxsubc^d;
2051 dst_slit=abs(xcent(1)-xslit)
2052 if (dst_slit<=half*slit_width)
then
2053 ixs=floor((xcent(2)-(xs(1)-half*dxsg))/dxsg)+1
2055 ixsmax=min(ixs+3,numxs)
2056 iwl=floor((wlc-(wl(1)-half*dwlg))/dwlg)+1
2057 nwl=3*ceiling(wlwd/dwlg+1)
2058 iwlmin=max(1,iwl-nwl)
2059 iwlmax=min(iwl+nwl,numwl)
2061 do iwl=iwlmin,iwlmax
2062 do ixs=ixsmin,ixsmax
2063 xerfmin1=(wl(iwl)-half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
2064 xerfmax1=(wl(iwl)+half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
2065 xerfmin2=(xs(ixs)-half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
2066 xerfmax2=(xs(ixs)+half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
2067 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
2068 spectra(iwl,ixs)=spectra(iwl,ixs)+fluxsubc*factor
2078 deallocate(flux,v,pth,te)
2086 integer,
intent(in) :: qunit
2088 character(20) :: datatype
2091 character (30) :: ion
2092 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
2093 double precision :: t0,t1
2096 datatype=
'image_euv'
2101 print *,
'###################################################'
2102 print *,
'Systhesizing EUV image'
2103 write(*,
'(a,f8.3,a)')
' Wavelength: ',linecent,
' Angstrom'
2104 print *,
'Unit of EUV flux: DN s^-1 pixel^-1'
2109 call mpistop(
'EUV dat-resolution needs Cartesian or spherical native rays')
2111 print *,
'Data-resolution image requested; native output pixel sizes are reported below.'
2113 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2115 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2129 call mpistop(
'ERROR: Wrong LOS for synthesizing emission!')
2133 write(*,
'(a,f7.1,a,f7.1,a,f5.1,a,f5.1,a)')
' Pixel: ',spacersl*725.0,
' km x ',spacersl*725.0,
' km (', &
2134 spacersl,
' arcsec x ', spacersl,
' arcsec)'
2136 print *,
'Unit of length: arcsec (~725 km)'
2139 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2141 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2147 '] of the simulation box is located at [X=0,Y=0] of the image'
2150 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'
2153 call mpistop(
"EUV synthesis: this coordinate is not supported!")
2158 if (
mype==0) print *,
'time comsuming: ',t1-t0,
' s'
2159 if (
mype==0) print *,
'###################################################'
2166 integer,
intent(in) :: qunit
2168 character(20) :: datatype
2169 double precision :: RHESSI_rsl
2170 double precision :: t0,t1
2173 datatype=
'image_sxr'
2178 print *,
'###################################################'
2179 print *,
'Systhesizing SXR image (observed at 1 AU).'
2186 print *,
'Unit of SXR flux: photons cm^-2 s^-1 pixel^-1'
2187 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, &
2188 ' Mm (', rhessi_rsl,
' arcsec x ', rhessi_rsl,
' arcsec)'
2190 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2192 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2202 call mpistop(
'ERROR: Wrong LOS for synthesizing emission!')
2206 print *,
'Unit of SXR flux: photons cm^-2 s^-1 pixel^-1'
2207 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, &
2208 ' Mm (', rhessi_rsl,
' arcsec x ', rhessi_rsl,
' arcsec)'
2210 print *,
'Unit of length: arcsec (~725 km)'
2213 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2215 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2221 '] of the simulation box is located at [X=0,Y=0] of the image'
2224 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'
2227 call mpistop(
"SXR synthesis: this coordinate is not supported!")
2232 if (
mype==0) print *,
'time comsuming:',t1-t0
2233 if (
mype==0) print *,
'###################################################'
2240 integer,
intent(in) :: qunit
2242 character(20) :: datatype
2243 double precision :: LASCO_rsl
2245 if (
mype==0) print *,
'###################################################'
2249 if (
mype==0) print *,
'Systhesizing white light image (observed by LASCO/C1).'
2252 if (
mype==0) print *,
'Systhesizing white light image (observed by LASCO/C2).'
2255 if (
mype==0) print *,
'Systhesizing white light image (observed by LASCO/C3).'
2257 call mpistop(
'Whitelight synthesis: instrument is not supported!')
2260 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 (', &
2261 lasco_rsl,
' arcsec x ', lasco_rsl,
' arcsec) '
2262 if (
mype==0) print *,
'Unit of white light flux: average Sun brightness'
2264 datatype=
'image_whitelight'
2269 print *,
'Unit of length: arcsec (~725 km)'
2272 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2274 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2280 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'
2283 call mpistop(
"Whitelight synthesis: this coordinate is not supported!")
2286 if (
mype==0) print *,
'###################################################'
2291 EUV,Dpl,nOut1,nOut2,xOut1,xOut2,&
2292 dxOut1,dxOut2,wOut,numWOut,Tau,EUVthin)
2296 integer,
intent(in) :: nSrc1,nSrc2
2297 double precision,
intent(in) :: xSrc1(nSrc1),xSrc2(nSrc2)
2298 double precision,
intent(in) :: dxSrc1(nSrc1),dxSrc2(nSrc2)
2299 double precision,
intent(in) :: EUV(nSrc1,nSrc2),Dpl(nSrc1,nSrc2)
2300 integer,
intent(out) :: nOut1,nOut2,numWOut
2301 double precision,
allocatable,
intent(out) :: xOut1(:),xOut2(:),dxOut1(:),dxOut2(:)
2302 double precision,
allocatable,
intent(out) :: wOut(:,:,:)
2303 double precision,
intent(in),
optional :: Tau(nSrc1,nSrc2),EUVthin(nSrc1,nSrc2)
2305 integer :: mass,ixS1,ixS2,ixP1,ixP2,ixC1,ixC2,iw
2306 integer :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
2307 character(30) :: ion
2308 double precision :: logTe,lineCent,spaceRsl,wlRsl,sigma_PSF,wslit
2309 double precision :: arcsec,dxInst,xMin1,xMax1,xMin2,xMax2,xCent1,xCent2
2310 double precision :: sigma0,xerfmin1,xerfmax1,xerfmin2,xerfmax2
2311 double precision :: factor,weightSum,weightNorm,thinVal,tauVal
2312 double precision,
allocatable :: dplNum(:,:),thinOut(:,:),tauOut(:,:),tauWeight(:,:)
2320 dxinst=spacersl*arcsec
2321 if (dxinst<=
zero)
call mpistop(
"instrument_postprocess has non-positive pixel size")
2323 xmin1=minval(xsrc1-
half*dxsrc1)
2324 xmax1=maxval(xsrc1+
half*dxsrc1)
2325 xmin2=minval(xsrc2-
half*dxsrc2)
2326 xmax2=maxval(xsrc2+
half*dxsrc2)
2327 xcent1=
half*(xmin1+xmax1)
2328 xcent2=
half*(xmin2+xmax2)
2329 nout1=16*max(1,ceiling((xmax1-xmin1)/(16.d0*dxinst)))
2330 nout2=16*max(1,ceiling((xmax2-xmin2)/(16.d0*dxinst)))
2331 xmin1=xcent1-
half*dble(nout1)*dxinst
2332 xmin2=xcent2-
half*dble(nout2)*dxinst
2334 allocate(xout1(nout1),xout2(nout2),dxout1(nout1),dxout2(nout2))
2336 xout1(ixp1)=xmin1+dxinst*(dble(ixp1)-
half)
2340 xout2(ixp2)=xmin2+dxinst*(dble(ixp2)-
half)
2345 if (
present(tau) .and.
output_tau) numwout=numwout+1
2347 allocate(wout(nout1,nout2,numwout),dplnum(nout1,nout2))
2350 if (
present(euvthin))
then
2351 allocate(thinout(nout1,nout2))
2355 allocate(tauout(nout1,nout2),tauweight(nout1,nout2))
2360 sigma0=sigma_psf*dxinst
2365 if (
present(euvthin)) thinval=euvthin(ixs1,ixs2)
2366 if (
present(tau)) tauval=tau(ixs1,ixs2)
2370 ixc1=floor((xsrc1(ixs1)-(xout1(1)-
half*dxinst))/dxinst)+1
2371 ixc2=floor((xsrc2(ixs2)-(xout2(1)-
half*dxinst))/dxinst)+1
2372 ixpmin1=max(1,ixc1-3)
2373 ixpmax1=min(nout1,ixc1+3)
2374 ixpmin2=max(1,ixc2-3)
2375 ixpmax2=min(nout2,ixc2+3)
2378 do ixp1=ixpmin1,ixpmax1
2379 do ixp2=ixpmin2,ixpmax2
2380 xerfmin1=((xout1(ixp1)-
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2381 xerfmax1=((xout1(ixp1)+
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2382 xerfmin2=((xout2(ixp2)-
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2383 xerfmax2=((xout2(ixp2)+
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2384 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
2385 weightsum=weightsum+factor
2388 if (weightsum<=
zero) cycle
2390 do ixp1=ixpmin1,ixpmax1
2391 do ixp2=ixpmin2,ixpmax2
2392 xerfmin1=((xout1(ixp1)-
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2393 xerfmax1=((xout1(ixp1)+
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2394 xerfmin2=((xout2(ixp2)-
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2395 xerfmax2=((xout2(ixp2)+
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2396 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
2397 weightnorm=factor/weightsum
2398 wout(ixp1,ixp2,1)=wout(ixp1,ixp2,1)+euv(ixs1,ixs2)*weightnorm
2399 dplnum(ixp1,ixp2)=dplnum(ixp1,ixp2)+euv(ixs1,ixs2)*dpl(ixs1,ixs2)*weightnorm
2400 if (
present(euvthin)) thinout(ixp1,ixp2)=thinout(ixp1,ixp2)+thinval*weightnorm
2402 tauout(ixp1,ixp2)=tauout(ixp1,ixp2)+tauval*weightnorm
2403 tauweight(ixp1,ixp2)=tauweight(ixp1,ixp2)+weightnorm
2413 wout(ixp1,ixp2,2)=dplnum(ixp1,ixp2)/wout(ixp1,ixp2,1)
2415 wout(ixp1,ixp2,2)=
zero
2424 if (tauweight(ixp1,ixp2)>
zero)
then
2425 wout(ixp1,ixp2,iw)=tauout(ixp1,ixp2)/tauweight(ixp1,ixp2)
2427 wout(ixp1,ixp2,iw)=
zero
2437 wout(ixp1,ixp2,iw)=min(
one,max(
zero,(thinout(ixp1,ixp2)-wout(ixp1,ixp2,1))/thinout(ixp1,ixp2)))
2439 wout(ixp1,ixp2,iw)=
zero
2446 write(*,
'(a,2(i8,1x),a,2(i8,1x),a,1pe12.5)') &
2447 ' instrument_postprocess EUV grid src/out: ',nsrc1,nsrc2,
' -> ',nout1,nout2,
' dx=',dxinst
2451 if (
allocated(thinout))
deallocate(thinout)
2452 if (
allocated(tauout))
deallocate(tauout,tauweight)
2456 Bright,nOut1,nOut2,xOut1,xOut2,&
2457 dxOut1,dxOut2,wOut,numWOut,Tau,BrightThin)
2460 integer,
intent(in) :: nSrc1,nSrc2
2461 double precision,
intent(in) :: xSrc1(nSrc1),xSrc2(nSrc2)
2462 double precision,
intent(in) :: dxSrc1(nSrc1),dxSrc2(nSrc2)
2463 double precision,
intent(in) :: Bright(nSrc1,nSrc2)
2464 integer,
intent(out) :: nOut1,nOut2,numWOut
2465 double precision,
allocatable,
intent(out) :: xOut1(:),xOut2(:),dxOut1(:),dxOut2(:)
2466 double precision,
allocatable,
intent(out) :: wOut(:,:,:)
2467 double precision,
intent(in),
optional :: Tau(nSrc1,nSrc2),BrightThin(nSrc1,nSrc2)
2469 integer :: ixS1,ixS2,ixP1,ixP2,ixC1,ixC2,iw,nStencil
2470 integer :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
2471 double precision :: arcsec,beamPixel,beamSigma,xMin1,xMax1,xMin2,xMax2,xCent1,xCent2
2472 double precision :: distance1,distance2,weight,cellArea,thinVal,tauVal
2473 double precision,
allocatable :: norm(:,:),thinOut(:,:),tauOut(:,:),tauNorm(:,:)
2486 if (beamsigma<=zero .or. beampixel<=zero)
then
2487 call mpistop(
"radio beam postprocess has non-positive beam or pixel size")
2490 xmin1=minval(xsrc1-half*dxsrc1)
2491 xmax1=maxval(xsrc1+half*dxsrc1)
2492 xmin2=minval(xsrc2-half*dxsrc2)
2493 xmax2=maxval(xsrc2+half*dxsrc2)
2494 xcent1=half*(xmin1+xmax1)
2495 xcent2=half*(xmin2+xmax2)
2496 nout1=16*max(1,ceiling((xmax1-xmin1)/(16.d0*beampixel)))
2497 nout2=16*max(1,ceiling((xmax2-xmin2)/(16.d0*beampixel)))
2498 xmin1=xcent1-half*dble(nout1)*beampixel
2499 xmin2=xcent2-half*dble(nout2)*beampixel
2501 allocate(xout1(nout1),xout2(nout2),dxout1(nout1),dxout2(nout2))
2503 xout1(ixp1)=xmin1+beampixel*(dble(ixp1)-half)
2504 dxout1(ixp1)=beampixel
2507 xout2(ixp2)=xmin2+beampixel*(dble(ixp2)-half)
2508 dxout2(ixp2)=beampixel
2512 if (
present(tau) .and.
output_tau) numwout=numwout+1
2514 allocate(wout(nout1,nout2,numwout),norm(nout1,nout2))
2517 if (
present(brightthin))
then
2518 allocate(thinout(nout1,nout2))
2522 allocate(tauout(nout1,nout2),taunorm(nout1,nout2))
2527 nstencil=max(3,ceiling(4.d0*beamsigma/beampixel)+1)
2532 if (
present(brightthin)) thinval=brightthin(ixs1,ixs2)
2533 if (
present(tau)) tauval=tau(ixs1,ixs2)
2534 if (abs(bright(ixs1,ixs2))<=smalldouble .and. abs(thinval)<=smalldouble .and. &
2535 abs(tauval)<=smalldouble) cycle
2537 ixc1=floor((xsrc1(ixs1)-(xout1(1)-half*beampixel))/beampixel)+1
2538 ixc2=floor((xsrc2(ixs2)-(xout2(1)-half*beampixel))/beampixel)+1
2539 ixpmin1=max(1,ixc1-nstencil)
2540 ixpmax1=min(nout1,ixc1+nstencil)
2541 ixpmin2=max(1,ixc2-nstencil)
2542 ixpmax2=min(nout2,ixc2+nstencil)
2543 cellarea=max(smalldouble,dxsrc1(ixs1)*dxsrc2(ixs2))
2545 do ixp1=ixpmin1,ixpmax1
2546 distance1=xout1(ixp1)-xsrc1(ixs1)
2547 do ixp2=ixpmin2,ixpmax2
2548 distance2=xout2(ixp2)-xsrc2(ixs2)
2549 weight=
exp_clamped(-half*(distance1**2+distance2**2)/beamsigma**2)*cellarea
2550 wout(ixp1,ixp2,1)=wout(ixp1,ixp2,1)+bright(ixs1,ixs2)*weight
2551 norm(ixp1,ixp2)=norm(ixp1,ixp2)+weight
2552 if (
present(brightthin)) thinout(ixp1,ixp2)=thinout(ixp1,ixp2)+thinval*weight
2554 tauout(ixp1,ixp2)=tauout(ixp1,ixp2)+tauval*weight
2555 taunorm(ixp1,ixp2)=taunorm(ixp1,ixp2)+weight
2564 if (norm(ixp1,ixp2)>zero)
then
2565 wout(ixp1,ixp2,1)=wout(ixp1,ixp2,1)/norm(ixp1,ixp2)
2566 if (
present(brightthin)) thinout(ixp1,ixp2)=thinout(ixp1,ixp2)/norm(ixp1,ixp2)
2568 wout(ixp1,ixp2,1)=zero
2569 if (
present(brightthin)) thinout(ixp1,ixp2)=zero
2579 if (taunorm(ixp1,ixp2)>zero)
then
2580 wout(ixp1,ixp2,iw)=tauout(ixp1,ixp2)/taunorm(ixp1,ixp2)
2582 wout(ixp1,ixp2,iw)=zero
2591 if (thinout(ixp1,ixp2)>smalldouble)
then
2592 wout(ixp1,ixp2,iw)=min(one,max(zero,(thinout(ixp1,ixp2)-wout(ixp1,ixp2,1))/thinout(ixp1,ixp2)))
2594 wout(ixp1,ixp2,iw)=zero
2601 write(*,
'(a,2(i8,1x),a,2(i8,1x),a,2(1pe12.5,1x))') &
2602 ' radio_beam_postprocess grid src/out: ',nsrc1,nsrc2,
' -> ',nout1,nout2,&
2607 if (
allocated(thinout))
deallocate(thinout)
2608 if (
allocated(tauout))
deallocate(tauout,taunorm)
2617 integer,
intent(in) :: qunit
2618 character(20),
intent(in) :: datatype
2621 double precision :: dx^D
2622 integer :: numX^D,ix^D
2623 double precision,
allocatable :: EUV(:,:),EUVs(:,:),Dpl(:,:),Dpls(:,:)
2624 double precision,
allocatable :: EUVthin(:,:),Tau(:,:)
2625 double precision,
allocatable :: SXR(:,:),SXRs(:,:),wI(:,:,:)
2626 double precision,
allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:),dxIi
2627 integer :: numXI1,numXI2,numSI,numWI,iw
2628 double precision :: xI^L
2629 integer :: iigrid,igrid,i,j
2630 double precision,
allocatable :: xIF1(:),xIF2(:),dxIF1(:),dxIF2(:)
2631 double precision,
allocatable :: xIP1(:),xIP2(:),dxIP1(:),dxIP2(:),wIP(:,:,:)
2632 integer :: nXIF1,nXIF2
2633 integer :: nXIP1,nXIP2,numWIP
2634 double precision :: xIF^L
2635 double precision :: vec_cor(1:3),xI_cor(1:2),dxDDA,xIcent1,xIcent2
2637 double precision :: unitv,arcsec,RHESSI_rsl,length_to_km
2638 integer :: strtype^D,nstrb^D,nbb^D,nuni^D,nstr^D,bnx^D
2639 double precision :: qs^D,dxfirst^D,dxmid^D,lenstr^D
2640 logical :: has_doppler_output,has_thick_output
2653 xicent1=
half*(xifmin1+xifmax1)
2654 xicent2=
half*(xifmin2+xifmax2)
2655 nxif1=max(1,ceiling((xifmax1-xifmin1)/dxdda))
2656 nxif2=max(1,ceiling((xifmax2-xifmin2)/dxdda))
2657 xifmin1=xicent1-
half*dble(nxif1)*dxdda
2658 xifmax1=xicent1+
half*dble(nxif1)*dxdda
2659 xifmin2=xicent2-
half*dble(nxif2)*dxdda
2660 xifmax2=xicent2+
half*dble(nxif2)*dxdda
2671 if (
mype==0)
write(*,
'(a,a,a,1pe12.5,a,2(i8,1x))') &
2673 ' image-plane dx=',dxdda,
' n=',nxif1,nxif2
2676 if (ix1==1) vec_cor(1)=xprobmin1
2677 if (ix1==2) vec_cor(1)=xprobmax1
2679 if (ix2==1) vec_cor(2)=xprobmin2
2680 if (ix2==2) vec_cor(2)=xprobmax2
2682 if (ix3==1) vec_cor(3)=xprobmin3
2683 if (ix3==2) vec_cor(3)=xprobmax3
2685 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
2691 xifmin1=min(xifmin1,xi_cor(1))
2692 xifmax1=max(xifmax1,xi_cor(1))
2693 xifmin2=min(xifmin2,xi_cor(2))
2694 xifmax2=max(xifmax2,xi_cor(2))
2700 xicent1=
half*(xifmin1+xifmax1)
2701 xicent2=
half*(xifmin2+xifmax2)
2702 nxif1=max(1,ceiling((xifmax1-xifmin1)/dxdda))
2703 nxif2=max(1,ceiling((xifmax2-xifmin2)/dxdda))
2704 xifmin1=xicent1-
half*dble(nxif1)*dxdda
2705 xifmax1=xicent1+
half*dble(nxif1)*dxdda
2706 xifmin2=xicent2-
half*dble(nxif2)*dxdda
2707 xifmax2=xicent2+
half*dble(nxif2)*dxdda
2718 if (
mype==0)
write(*,
'(a,a,a,1pe12.5,a,2(i8,1x))') &
2720 ' image-plane dx=',dxdda,
' n=',nxif1,nxif2
2738 if (
mype==0)
write(*,
'(a)')
' LOS vector: [-1.00 0.00 0.00]'
2739 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [ 0.00 1.00 0.00]'
2740 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 0.00 1.00]'
2758 if (
mype==0)
write(*,
'(a)')
' LOS vector: [ 0.00 -1.00 0.00]'
2759 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [-1.00 0.00 0.00]'
2760 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 0.00 1.00]'
2778 if (
mype==0)
write(*,
'(a)')
' LOS vector: [ 0.00 0.00 -1.00]'
2779 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [ 1.00 0.00 0.00]'
2780 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 1.00 0.00]'
2782 allocate(xif1(nxif1),xif2(nxif2),dxif1(nxif1),dxif2(nxif2))
2785 select case(strtype1)
2787 dxif1(:)=(xifmax1-xifmin1)/nxif1
2789 xif1(ix1)=xifmin1+dxif1(ix1)*(ix1-
half)
2793 dxfirst1=(xifmax1-xifmin1)*(
one-qs1)/(
one-qs1**nxif1)
2796 dxif1(ix1)=dxfirst1*qs1**(ix1-1)
2797 xif1(ix1)=dxif1(1)/(
one-qs1)*(
one-qs1**(ix1-1))+
half*dxif1(ix1)
2802 nuni1=nbb1-nstrb1*bnx1
2803 lenstr1=(xifmax1-xifmin1)/(2.d0+nuni1*(
one-qs1)/(
one-qs1**nstr1))
2804 dxfirst1=(xifmax1-xifmin1)/(dble(nuni1)+2.d0/(
one-qs1)*(
one-qs1**nstr1))
2810 dxfirst1=lenstr1*(
one-qs1)/(
one-qs1**nstr1)
2813 if(nuni1 .gt. 0)
then
2814 do ix1=nstr1+1,nstr1+nuni1
2816 xif1(ix1)=lenstr1+(dble(ix1)-0.5d0-nstr1)*dxif1(ix1)+xifmin1
2821 dxif1(ix1)=dxfirst1*qs1**(nstr1-ix1)
2822 xif1(ix1)=xifmin1+lenstr1-dxif1(ix1)*
half-dxfirst1*(
one-qs1**(nstr1-ix1))/(
one-qs1)
2825 do ix1=nstr1+nuni1+1,nxif1
2826 dxif1(ix1)=dxfirst1*qs1**(ix1-nstr1-nuni1-1)
2827 xif1(ix1)=xifmax1-lenstr1+dxif1(ix1)*
half+dxfirst1*(
one-qs1**(ix1-nstr1-nuni1-1))/(
one-qs1)
2830 call mpistop(
"unknown stretch type")
2833 select case(strtype2)
2835 dxif2(:)=(xifmax2-xifmin2)/nxif2
2837 xif2(ix2)=xifmin2+dxif2(ix2)*(ix2-
half)
2841 dxfirst2=(xifmax2-xifmin2)*(
one-qs2)/(
one-qs2**nxif2)
2844 dxif2(ix2)=dxfirst2*qs2**(ix2-1)
2845 xif2(ix2)=dxif2(1)/(
one-qs2)*(
one-qs2**(ix2-1))+
half*dxif2(ix2)
2850 nuni2=nbb2-nstrb2*bnx2
2851 lenstr2=(xifmax2-xifmin2)/(2.d0+nuni2*(
one-qs2)/(
one-qs2**nstr2))
2852 dxfirst2=(xifmax2-xifmin2)/(dble(nuni2)+2.d0/(
one-qs2)*(
one-qs2**nstr2))
2858 dxfirst2=lenstr2*(
one-qs2)/(
one-qs2**nstr2)
2861 if(nuni2 .gt. 0)
then
2862 do ix2=nstr2+1,nstr2+nuni2
2864 xif2(ix2)=lenstr2+(dble(ix2)-0.5d0-nstr2)*dxif2(ix2)+xifmin2
2869 dxif2(ix2)=dxfirst2*qs2**(nstr2-ix2)
2870 xif2(ix2)=xifmin2+lenstr2-dxif2(ix2)*
half-dxfirst2*(
one-qs2**(nstr2-ix2))/(
one-qs2)
2873 do ix2=nstr2+nuni2+1,nxif2
2874 dxif2(ix2)=dxfirst2*qs2**(ix2-nstr2-nuni2-1)
2875 xif2(ix2)=xifmax2-lenstr2+dxif2(ix2)*
half+dxfirst2*(
one-qs2**(ix2-nstr2-nuni2-1))/(
one-qs2)
2878 call mpistop(
"unknown stretch type")
2881 if (
mype==0 .and. datatype==
'image_euv')
then
2889 write(*,
'(a,i8,a,i8)')
' Native data-resolution image grid: ',nxif1,
' x ',nxif2
2890 write(*,
'(a,f10.3,a,f10.3,a,f8.3,a,f8.3,a)') &
2891 ' Native xI1 pixel-size range: ',minval(dxif1)*length_to_km,
'--', &
2892 maxval(dxif1)*length_to_km,
' km (',minval(dxif1)/arcsec,
'--',maxval(dxif1)/arcsec,
' arcsec)'
2893 write(*,
'(a,f10.3,a,f10.3,a,f8.3,a,f8.3,a)') &
2894 ' Native xI2 pixel-size range: ',minval(dxif2)*length_to_km,
'--', &
2895 maxval(dxif2)*length_to_km,
' km (',minval(dxif2)/arcsec,
'--',maxval(dxif2)/arcsec,
' arcsec)'
2899 if (datatype==
'image_euv')
then
2908 allocate(wi(nxif1,nxif2,numwi))
2909 allocate(euv(nxif1,nxif2),dpl(nxif1,nxif2))
2911 allocate(euvthin(nxif1,nxif2),tau(nxif1,nxif2))
2920 if (has_doppler_output)
then
2925 allocate(euvs(nxif1,nxif2),dpls(nxif1,nxif2))
2934 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision, &
2939 do iigrid=1,igridstail; igrid=igrids(iigrid);
2943 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision, &
2945 call mpi_allreduce(dpls,dpl,numsi,mpi_double_precision, &
2948 if (has_doppler_output)
then
2952 deallocate(euvs,dpls)
2954 if (has_thick_output)
then
2955 if (has_doppler_output)
then
2957 has_thick_output,dpl=dpl,tau=tau,euvthin=euvthin)
2960 has_thick_output,tau=tau,euvthin=euvthin)
2962 else if (has_doppler_output)
then
2964 has_thick_output,dpl=dpl)
2974 euv,nxip1,nxip2,xip1,xip2,&
2975 dxip1,dxip2,wip,numwip,tau=tau,brightthin=euvthin)
2978 euv,nxip1,nxip2,xip1,xip2,&
2979 dxip1,dxip2,wip,numwip)
2983 euv,dpl,nxip1,nxip2,xip1,xip2,&
2984 dxip1,dxip2,wip,numwip,tau=tau,euvthin=euvthin)
2987 euv,dpl,nxip1,nxip2,xip1,xip2,&
2988 dxip1,dxip2,wip,numwip)
2990 call output_data(qunit,xip1,xip2,dxip1,dxip2,wip,nxip1,nxip2,numwip,datatype)
2991 deallocate(xip1,xip2,dxip1,dxip2,wip)
2993 call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
2996 deallocate(wi,euv,dpl,euvthin,tau)
2998 deallocate(wi,euv,dpl)
3003 if (datatype==
'image_sxr')
then
3011 allocate(wi(nxif1,nxif2,numwi))
3012 allocate(sxrs(nxif1,nxif2),sxr(nxif1,nxif2))
3015 do iigrid=1,igridstail; igrid=igrids(iigrid);
3019 call mpi_allreduce(sxrs,sxr,numsi,mpi_double_precision, &
3022 sxr=sxr*(rhessi_rsl*arcsec)**2
3030 call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
3031 deallocate(wi,sxr,sxrs)
3034 deallocate(xif1,xif2,dxif1,dxif2)
3041 integer,
intent(in) :: igrid,nXIF1,nXIF2
3042 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3043 double precision,
intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
3045 double precision,
intent(out) :: SXR(nXIF1,nXIF2)
3047 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
3048 double precision :: xb^L,xd^D
3049 double precision,
allocatable :: flux(:^D&),opacity(:^D&)
3050 double precision,
allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
3051 double precision,
allocatable :: SXRg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
3052 integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
3053 double precision :: SXRt,xc^L,xg^L,r2,area_1AU
3054 integer :: ixP^L,ixP^D
3055 integer :: direction_LOS
3065 ^d&ixomin^d=ixmlo^d\
3066 ^d&ixomax^d=ixmhi^d\
3067 ^d&iximin^d=
ixglo^d\
3068 ^d&iximax^d=
ixghi^d\
3072 allocate(flux(ixi^s))
3073 allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
3074 dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
3075 dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
3076 dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
3081 levelg=ps(igrid)%level
3085 select case(direction_los)
3096 allocate(sxrg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
3102 select case(direction_los)
3104 do ix2=ixomin2,ixomax2
3105 ixgmin1=(ix2-1)*rft+1
3107 do ix3=ixomin3,ixomax3
3108 ixgmin2=(ix3-1)*rft+1
3111 do ix1=ixomin1,ixomax1
3114 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
3118 do ix3=ixomin3,ixomax3
3119 ixgmin1=(ix3-1)*rft+1
3121 do ix1=ixomin1,ixomax1
3122 ixgmin2=(ix1-1)*rft+1
3125 do ix2=ixomin2,ixomax2
3128 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
3132 do ix1=ixomin1,ixomax1
3133 ixgmin1=(ix1-1)*rft+1
3135 do ix2=ixomin2,ixomax2
3136 ixgmin2=(ix2-1)*rft+1
3139 do ix3=ixomin3,ixomax3
3142 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
3152 select case(direction_los)
3154 ixgmin1=(ixomin2-1)*rft+1
3156 ixgmin2=(ixomin3-1)*rft+1
3159 ixgmin1=(ixomin3-1)*rft+1
3161 ixgmin2=(ixomin1-1)*rft+1
3164 ixgmin1=(ixomin1-1)*rft+1
3166 ixgmin2=(ixomin2-1)*rft+1
3170 select case(direction_los)
3172 ixpmin1=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3173 ixpmax1=
node(pig2_,igrid)*rft*block_nx2
3174 ixpmin2=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3175 ixpmax2=
node(pig3_,igrid)*rft*block_nx3
3177 ixpmin1=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3178 ixpmax1=
node(pig3_,igrid)*rft*block_nx3
3179 ixpmin2=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3180 ixpmax2=
node(pig1_,igrid)*rft*block_nx1
3182 ixpmin1=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3183 ixpmax1=
node(pig1_,igrid)*rft*block_nx1
3184 ixpmin2=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3185 ixpmax2=
node(pig2_,igrid)*rft*block_nx2
3187 xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
3188 xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
3189 dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
3190 dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
3191 sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
3192 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
3194 deallocate(flux,dxb1,dxb2,dxb3,sxrg,xg1,xg2,dxg1,dxg2)
3201 integer,
intent(in) :: igrid,nXIF1,nXIF2
3202 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3203 double precision,
intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
3205 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
3207 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
3208 double precision :: xb^L,xd^D
3209 double precision,
allocatable :: flux(:^D&),v(:^D&),rho(:^D&),opacity(:^D&)
3210 double precision,
allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
3211 double precision,
allocatable :: EUVg(:,:),Fvg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
3212 integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
3213 double precision :: EUVt,Fvt,xc^L,xg^L,r2
3214 integer :: ixP^L,ixP^D
3215 integer :: direction_LOS
3225 ^d&ixomin^d=ixmlo^d\
3226 ^d&ixomax^d=ixmhi^d\
3227 ^d&iximin^d=
ixglo^d\
3228 ^d&iximax^d=
ixghi^d\
3232 allocate(flux(ixi^s),v(ixi^s),rho(ixi^s),opacity(ixi^s))
3233 allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
3234 dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
3235 dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
3236 dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
3249 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
3250 v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
3255 levelg=ps(igrid)%level
3259 select case(direction_los)
3270 allocate(euvg(nxg1,nxg2),fvg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
3277 select case(direction_los)
3279 do ix2=ixomin2,ixomax2
3280 ixgmin1=(ix2-1)*rft+1
3282 do ix3=ixomin3,ixomax3
3283 ixgmin2=(ix3-1)*rft+1
3287 do ix1=ixomin1,ixomax1
3291 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
3292 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
3296 do ix3=ixomin3,ixomax3
3297 ixgmin1=(ix3-1)*rft+1
3299 do ix1=ixomin1,ixomax1
3300 ixgmin2=(ix1-1)*rft+1
3304 do ix2=ixomin2,ixomax2
3308 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
3309 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
3313 do ix1=ixomin1,ixomax1
3314 ixgmin1=(ix1-1)*rft+1
3316 do ix2=ixomin2,ixomax2
3317 ixgmin2=(ix2-1)*rft+1
3321 do ix3=ixomin3,ixomax3
3325 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
3326 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
3337 select case(direction_los)
3339 ixgmin1=(ixomin2-1)*rft+1
3341 ixgmin2=(ixomin3-1)*rft+1
3344 ixgmin1=(ixomin3-1)*rft+1
3346 ixgmin2=(ixomin1-1)*rft+1
3349 ixgmin1=(ixomin1-1)*rft+1
3351 ixgmin2=(ixomin2-1)*rft+1
3355 select case(direction_los)
3357 ixpmin1=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3358 ixpmax1=
node(pig2_,igrid)*rft*block_nx2
3359 ixpmin2=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3360 ixpmax2=
node(pig3_,igrid)*rft*block_nx3
3362 ixpmin1=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3363 ixpmax1=
node(pig3_,igrid)*rft*block_nx3
3364 ixpmin2=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3365 ixpmax2=
node(pig1_,igrid)*rft*block_nx1
3367 ixpmin1=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3368 ixpmax1=
node(pig1_,igrid)*rft*block_nx1
3369 ixpmin2=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3370 ixpmax2=
node(pig2_,igrid)*rft*block_nx2
3372 xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
3373 xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
3374 dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
3375 dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
3376 euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
3377 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
3378 dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
3379 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
3381 deallocate(flux,v,opacity,dxb1,dxb2,dxb3,euvg,fvg,xg1,xg2,dxg1,dxg2)
3390 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),box_min(1:3),box_max(1:3)
3391 logical,
intent(out) :: hit
3392 double precision,
intent(out) :: t_enter,t_exit
3395 double precision :: t1,t2,td
3401 if (abs(ray_dir(idir))<=smalldouble)
then
3402 if (ray_origin(idir)<box_min(idir) .or. ray_origin(idir)>box_max(idir))
then
3407 t1=(box_min(idir)-ray_origin(idir))/ray_dir(idir)
3408 t2=(box_max(idir)-ray_origin(idir))/ray_dir(idir)
3414 t_enter=max(t_enter,t1)
3415 t_exit=min(t_exit,t2)
3416 if (t_enter>=t_exit)
then
3425 integer,
intent(in) :: ixI^L, ixO^L
3426 double precision,
intent(in) :: x(ixI^S,1:ndim),dx(ixI^S,1:ndim)
3427 double precision,
allocatable,
intent(out) :: xface1(:),xface2(:),xface3(:)
3431 allocate(xface1(ixomin1:ixomax1+1),xface2(ixomin2:ixomax2+1),xface3(ixomin3:ixomax3+1))
3435 do ix1=ixomin1,ixomax1
3436 xface1(ix1)=x(ix^d,1)-half*dx(ix^d,1)
3438 xface1(ixomax1+1)=x(ixomax1,ixomin2,ixomin3,1)+half*dx(ixomax1,ixomin2,ixomin3,1)
3442 do ix2=ixomin2,ixomax2
3443 xface2(ix2)=x(ix^d,2)-half*dx(ix^d,2)
3445 xface2(ixomax2+1)=x(ixomin1,ixomax2,ixomin3,2)+half*dx(ixomin1,ixomax2,ixomin3,2)
3449 do ix3=ixomin3,ixomax3
3450 xface3(ix3)=x(ix^d,3)-half*dx(ix^d,3)
3452 xface3(ixomax3+1)=x(ixomin1,ixomin2,ixomax3,3)+half*dx(ixomin1,ixomin2,ixomax3,3)
3456 integer,
intent(in) :: imin,imax
3457 double precision,
intent(in) :: pos,faces(imin:imax+1)
3459 integer :: ilo,ihi,imid
3461 if (pos<=faces(imin))
then
3465 if (pos>=faces(imax+1))
then
3472 do while (ihi-ilo>1)
3474 if (pos>=faces(imid))
then
3480 idx=min(imax,max(imin,ilo))
3484 integer,
intent(in) :: imin,imax,idx
3485 double precision,
intent(in) :: ray_origin_axis,ray_dir_axis,faces(imin:imax+1)
3486 integer,
intent(out) :: step
3487 double precision,
intent(out) :: tMax
3489 if (ray_dir_axis>zero)
then
3491 tmax=(faces(idx+1)-ray_origin_axis)/ray_dir_axis
3492 else if (ray_dir_axis<zero)
then
3494 tmax=(faces(idx)-ray_origin_axis)/ray_dir_axis
3502 integer,
intent(in) :: imin,imax,step
3503 double precision,
intent(in) :: ray_origin_axis,ray_dir_axis,faces(imin:imax+1)
3504 integer,
intent(inout) :: idx
3505 double precision,
intent(inout) :: tMax
3506 logical,
intent(out) :: done
3510 if (idx<imin .or. idx>imax)
then
3515 tmax=(faces(idx+1)-ray_origin_axis)/ray_dir_axis
3516 else if (step<0)
then
3517 tmax=(faces(idx)-ray_origin_axis)/ray_dir_axis
3524 ray_origin,xface1,xface2,xface3,t_enter,t_exit,EUVp,Dplp)
3525 integer,
intent(in) :: ixI^L, ixO^L
3526 double precision,
intent(in) :: source(ixI^S),sourcev(ixI^S)
3527 double precision,
intent(in) :: ray_origin(1:3)
3528 double precision,
intent(in) :: xface1(ixOmin1:ixOmax1+1),xface2(ixOmin2:ixOmax2+1),&
3529 xface3(ixOmin3:ixOmax3+1)
3530 double precision,
intent(in) :: t_enter,t_exit
3531 double precision,
intent(inout) :: EUVp,Dplp
3533 integer :: ix^D,step(1:3)
3534 double precision :: pos(1:3),tMax(1:3),tNow,tNext,ds_cm,epsRay
3537 if (t_exit<=t_enter)
return
3538 epsray=max(1.d-12,1.d-10*abs(t_exit-t_enter))
3539 pos=ray_origin+(t_enter+epsray)*
vec_los
3550 tnext=min(t_exit,tmax(1),tmax(2),tmax(3))
3551 if (tnext>tnow)
then
3552 ds_cm=(tnext-tnow)*unit_length
3553 if (si_unit) ds_cm=ds_cm*1.d2
3554 euvp=euvp+source(ix^d)*ds_cm
3555 dplp=dplp+sourcev(ix^d)*ds_cm
3558 if (tnow>=t_exit-epsray)
exit
3560 if (tmax(1)<=tnow+epsray)
then
3564 if (tmax(2)<=tnow+epsray)
then
3568 if (tmax(3)<=tnow+epsray)
then
3576 double precision,
allocatable,
intent(inout) :: segments(:,:)
3577 integer,
intent(inout) :: nseg,capacity
3578 integer,
intent(in) :: pixel_id
3579 double precision,
intent(in) :: tseg,jds,kds,jvds
3581 double precision,
allocatable :: tmp(:,:)
3582 integer :: new_capacity
3584 if (capacity<=0)
then
3586 allocate(segments(5,capacity))
3587 else if (nseg>=capacity)
then
3588 new_capacity=2*capacity
3589 allocate(tmp(5,new_capacity))
3590 tmp(:,1:capacity)=segments(:,1:capacity)
3591 call move_alloc(tmp,segments)
3592 capacity=new_capacity
3596 segments(1,nseg)=dble(pixel_id)
3597 segments(2,nseg)=tseg
3598 segments(3,nseg)=jds
3599 segments(4,nseg)=kds
3600 segments(5,nseg)=jvds
3604 pixel_id,ray_origin,xface1,xface2,xface3,t_enter,t_exit,&
3605 segments,nseg,capacity)
3606 integer,
intent(in) :: ixI^L, ixO^L
3607 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S),sourcev(ixI^S)
3608 integer,
intent(in) :: pixel_id
3609 double precision,
intent(in) :: ray_origin(1:3)
3610 double precision,
intent(in) :: xface1(ixOmin1:ixOmax1+1),xface2(ixOmin2:ixOmax2+1),&
3611 xface3(ixOmin3:ixOmax3+1)
3612 double precision,
intent(in) :: t_enter,t_exit
3613 double precision,
allocatable,
intent(inout) :: segments(:,:)
3614 integer,
intent(inout) :: nseg,capacity
3616 integer :: ix^D,step(1:3)
3617 double precision :: pos(1:3),tMax(1:3),tNow,tNext,ds_cm,epsRay,tseg
3618 double precision :: jds,kds,jvds
3621 if (t_exit<=t_enter)
return
3622 epsray=max(1.d-12,1.d-10*abs(t_exit-t_enter))
3623 pos=ray_origin+(t_enter+epsray)*
vec_los
3634 tnext=min(t_exit,tmax(1),tmax(2),tmax(3))
3635 if (tnext>tnow)
then
3636 ds_cm=(tnext-tnow)*unit_length
3637 if (si_unit) ds_cm=ds_cm*1.d2
3638 jds=source(ix^d)*ds_cm
3639 kds=opacity(ix^d)*ds_cm
3640 jvds=sourcev(ix^d)*ds_cm
3641 if (jds/=zero .or. kds/=zero .or. jvds/=zero)
then
3642 tseg=half*(tnow+tnext)
3647 if (tnow>=t_exit-epsray)
exit
3649 if (tmax(1)<=tnow+epsray)
then
3653 if (tmax(2)<=tnow+epsray)
then
3657 if (tmax(3)<=tnow+epsray)
then
3665 double precision,
intent(in) :: segments(:,:)
3666 integer,
intent(inout) :: idx(:)
3667 integer,
intent(in) :: nidx
3678 double precision,
intent(in) :: segments(:,:)
3679 integer,
intent(inout) :: idx(:)
3680 integer,
intent(in) :: ilo,ihi
3687 do while (j>=ilo .and. segments(2,idx(j))>segments(2,key))
3696 double precision,
intent(in) :: segments(:,:)
3697 integer,
intent(inout) :: idx(:)
3698 integer,
intent(in) :: ilo,ihi
3701 double precision :: pivot
3703 if (ihi-ilo<=32)
then
3710 pivot=segments(2,idx((ilo+ihi)/2))
3712 do while (segments(2,idx(i))<pivot)
3715 do while (segments(2,idx(j))>pivot)
3733 integer,
intent(in) :: pixel_id
3735 owner=mod(pixel_id-1,npe)
3739 double precision,
intent(in) :: segments(:,:)
3740 integer,
intent(in) :: is,nvars
3746 if (segments(iv,is)/=segments(iv,is) .or. abs(segments(iv,is))>=1.d90)
then
3753 subroutine cart_dda_block_pixel_range(box_min,box_max,nXIF1,nXIF2,xIF1,xIF2,ixPmin1,ixPmax1,ixPmin2,ixPmax2,has_pixels)
3754 double precision,
intent(in) :: box_min(1:3),box_max(1:3)
3755 integer,
intent(in) :: nXIF1,nXIF2
3756 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3757 integer,
intent(out) :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
3758 logical,
intent(out) :: has_pixels
3761 double precision :: vec_cor(1:3),xI_cor(1:2)
3762 double precision :: xmin1,xmax1,xmin2,xmax2,dx1,dx2
3765 if (i1==1) vec_cor(1)=box_min(1)
3766 if (i1==2) vec_cor(1)=box_max(1)
3768 if (i2==1) vec_cor(2)=box_min(2)
3769 if (i2==2) vec_cor(2)=box_max(2)
3771 if (i3==1) vec_cor(3)=box_min(3)
3772 if (i3==2) vec_cor(3)=box_max(3)
3774 if (i1==1 .and. i2==1 .and. i3==1)
then
3780 xmin1=min(xmin1,xi_cor(1))
3781 xmax1=max(xmax1,xi_cor(1))
3782 xmin2=min(xmin2,xi_cor(2))
3783 xmax2=max(xmax2,xi_cor(2))
3790 dx1=abs(xif1(2)-xif1(1))
3792 dx1=max(one,abs(xmax1-xmin1))
3795 dx2=abs(xif2(2)-xif2(1))
3797 dx2=max(one,abs(xmax2-xmin2))
3800 ixpmin1=max(1,floor((xmin1-xif1(1))/dx1)+1-1)
3801 ixpmax1=min(nxif1,ceiling((xmax1-xif1(1))/dx1)+1+1)
3802 ixpmin2=max(1,floor((xmin2-xif2(1))/dx2)+1-1)
3803 ixpmax2=min(nxif2,ceiling((xmax2-xif2(1))/dx2)+1+1)
3804 has_pixels=ixpmin1<=ixpmax1 .and. ixpmin2<=ixpmax2
3810 integer,
intent(in) :: nXIF1,nXIF2
3811 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3813 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
3815 integer :: ixO^L,ixO^D,ixI^L,ix^D
3816 integer :: iigrid,igrid,ixP1,ixP2,numSI,ixPmin1,ixPmax1,ixPmin2,ixPmax2
3817 double precision :: box_min(1:3),box_max(1:3),ray_origin(1:3)
3818 double precision :: t_enter,t_exit,vlos
3819 double precision :: profile_local(2),profile_global(2)
3820 logical :: hit,has_pixels
3821 double precision,
allocatable :: source(:^D&),sourcev(:^D&),rho(:^D&),opacity(:^D&)
3822 double precision,
allocatable :: xface1(:),xface2(:),xface3(:)
3823 double precision,
allocatable :: EUVs(:,:),Dpls(:,:)
3825 allocate(euvs(nxif1,nxif2),dpls(nxif1,nxif2))
3830 do iigrid=1,igridstail; igrid=igrids(iigrid);
3831 ^d&ixomin^d=ixmlo^d\
3832 ^d&ixomax^d=ixmhi^d\
3833 ^d&iximin^d=
ixglo^d\
3834 ^d&iximax^d=
ixghi^d\
3836 box_min(1)=
rnode(rpxmin1_,igrid)
3837 box_min(2)=
rnode(rpxmin2_,igrid)
3838 box_min(3)=
rnode(rpxmin3_,igrid)
3839 box_max(1)=
rnode(rpxmax1_,igrid)
3840 box_max(2)=
rnode(rpxmax2_,igrid)
3841 box_max(3)=
rnode(rpxmax3_,igrid)
3844 ixpmin1,ixpmax1,ixpmin2,ixpmax2,has_pixels)
3845 if (.not. has_pixels)
then
3846 deallocate(xface1,xface2,xface3)
3850 allocate(source(ixi^s),sourcev(ixi^s),rho(ixi^s),opacity(ixi^s))
3860 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
3861 do ix1=ixomin1,ixomax1
3862 do ix2=ixomin2,ixomax2
3863 do ix3=ixomin3,ixomax3
3864 if (rho(ix^d)>smalldouble)
then
3865 vlos=(ps(igrid)%w(ix^d,iw_mom(1))*
vec_los(1)+&
3866 ps(igrid)%w(ix^d,iw_mom(2))*
vec_los(2)+&
3867 ps(igrid)%w(ix^d,iw_mom(3))*
vec_los(3))/rho(ix^d)
3868 sourcev(ix^d)=source(ix^d)*vlos
3874 deallocate(rho,opacity)
3876 do ixp1=ixpmin1,ixpmax1
3877 do ixp2=ixpmin2,ixpmax2
3879 profile_local(1)=profile_local(1)+one
3882 profile_local(2)=profile_local(2)+one
3884 ray_origin,xface1,xface2,xface3,t_enter,t_exit,euvs(ixp1,ixp2),dpls(ixp1,ixp2))
3889 deallocate(source,sourcev,xface1,xface2,xface3)
3893 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
3894 call mpi_allreduce(dpls,dpl,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
3895 call mpi_allreduce(profile_local,profile_global,2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
3897 write(*,
'(a,2(es12.5,1x))')
' cart_dda thin profile ray_tests ray_hits: ',profile_global
3899 deallocate(euvs,dpls)
3905 integer,
intent(in) :: nXIF1,nXIF2
3906 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3908 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
3909 double precision,
intent(out) :: Tau(nXIF1,nXIF2),EUVthin(nXIF1,nXIF2)
3911 integer,
parameter :: nSegVars=5
3912 integer :: ixO^L,ixO^D,ixI^L,ix^D
3913 integer :: iigrid,igrid,ixP1,ixP2,ipix,ipixStart,ipixEnd,nPixBatch,pixel_id
3914 integer :: nseg,capacity,totalCount,totalSeg,ipe,is,iseg,nidx,owner,isegDest,nsegBefore
3915 integer :: ixGlobal,iyGlobal,ixPmin1,ixPmax1,ixPmin2,ixPmax2,iFirst,iLast,iLocal
3916 integer :: nPixBatchTarget
3917 integer :: maxSegBatchTarget,maxSegCommTarget,maxNsegBatch,nPixTotal
3918 integer :: maxOwnerSegCount,maxOwnerSegCountLocal,segOffset,recvFill,totalRoundCount,totalRoundSeg
3919 integer,
allocatable :: sendCounts(:),recvCounts(:),sendDispls(:),recvDispls(:)
3920 integer,
allocatable :: roundSendCounts(:),roundRecvCounts(:)
3921 integer,
allocatable :: roundSendDispls(:),roundRecvDispls(:)
3922 integer,
allocatable :: ownerSegCounts(:),ownerOffsets(:),idx(:)
3923 integer,
allocatable :: bucketCounts(:),bucketOffsets(:),bucketFill(:)
3924 double precision :: ray_origin(1:3)
3925 double precision :: t_enter,t_exit,vlos,atten
3926 double precision :: profile_local(5),profile_global(5),profile_batch(5)
3927 logical :: hit,has_pixels,batchAccepted,batchReduced
3928 double precision,
allocatable :: rho(:^D&)
3929 double precision,
allocatable :: segments(:,:),segments_send(:,:),segments_recv(:,:)
3930 double precision,
allocatable :: segments_recv_round(:,:)
3931 double precision,
allocatable :: image_reduce(:,:)
3939 allocate(sendcounts(0:
npe-1),recvcounts(0:
npe-1),senddispls(0:
npe-1),recvdispls(0:
npe-1))
3940 allocate(roundsendcounts(0:
npe-1),roundrecvcounts(0:
npe-1))
3941 allocate(roundsenddispls(0:
npe-1),roundrecvdispls(0:
npe-1))
3942 allocate(ownersegcounts(0:
npe-1),owneroffsets(0:
npe-1))
3943 allocate(cache(igridstail))
3945 allocate(bucketcounts(npixbatchtarget),bucketoffsets(npixbatchtarget+1),&
3946 bucketfill(npixbatchtarget))
3948 do iigrid=1,igridstail; igrid=igrids(iigrid);
3949 ^d&ixomin^d=ixmlo^d\
3950 ^d&ixomax^d=ixmhi^d\
3951 ^d&iximin^d=
ixglo^d\
3952 ^d&iximax^d=
ixghi^d\
3954 cache(iigrid)%igrid=igrid
3955 allocate(cache(iigrid)%source(ixi^s),cache(iigrid)%opacity(ixi^s),&
3956 cache(iigrid)%sourcev(ixi^s),rho(ixi^s))
3957 cache(iigrid)%source=zero
3958 cache(iigrid)%opacity=zero
3959 cache(iigrid)%sourcev=zero
3962 cache(iigrid)%source,cache(iigrid)%opacity)
3964 call get_euv(
wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,cache(iigrid)%source)
3967 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
3968 do ix1=ixomin1,ixomax1
3969 do ix2=ixomin2,ixomax2
3970 do ix3=ixomin3,ixomax3
3971 if (rho(ix^d)>smalldouble)
then
3972 vlos=(ps(igrid)%w(ix^d,iw_mom(1))*
vec_los(1)+&
3973 ps(igrid)%w(ix^d,iw_mom(2))*
vec_los(2)+&
3974 ps(igrid)%w(ix^d,iw_mom(3))*
vec_los(3))/rho(ix^d)
3975 cache(iigrid)%sourcev(ix^d)=cache(iigrid)%source(ix^d)*vlos
3982 cache(iigrid)%box_min(1)=
rnode(rpxmin1_,igrid)
3983 cache(iigrid)%box_min(2)=
rnode(rpxmin2_,igrid)
3984 cache(iigrid)%box_min(3)=
rnode(rpxmin3_,igrid)
3985 cache(iigrid)%box_max(1)=
rnode(rpxmax1_,igrid)
3986 cache(iigrid)%box_max(2)=
rnode(rpxmax2_,igrid)
3987 cache(iigrid)%box_max(3)=
rnode(rpxmax3_,igrid)
3989 cache(iigrid)%xface1,cache(iigrid)%xface2,&
3990 cache(iigrid)%xface3)
3992 nxif1,nxif2,xif1,xif2,cache(iigrid)%ixPmin1,cache(iigrid)%ixPmax1,&
3993 cache(iigrid)%ixPmin2,cache(iigrid)%ixPmax2,cache(iigrid)%has_pixels)
3996 npixtotal=nxif1*nxif2
3998 do while (ipixstart<=npixtotal)
3999 ipixend=min(nxif1*nxif2,ipixstart+npixbatchtarget-1)
4000 npixbatch=ipixend-ipixstart+1
4001 batchaccepted=.false.
4002 batchreduced=.false.
4004 do while (.not. batchaccepted)
4009 do iigrid=1,igridstail; igrid=igrids(iigrid);
4010 ^d&ixomin^d=ixmlo^d\
4011 ^d&ixomax^d=ixmhi^d\
4012 ^d&iximin^d=
ixglo^d\
4013 ^d&iximax^d=
ixghi^d\
4015 ixpmin1=cache(iigrid)%ixPmin1
4016 ixpmax1=cache(iigrid)%ixPmax1
4017 ixpmin2=cache(iigrid)%ixPmin2
4018 ixpmax2=cache(iigrid)%ixPmax2
4019 has_pixels=cache(iigrid)%has_pixels
4020 if (.not. has_pixels) cycle
4022 do ixp2=ixpmin2,ixpmax2
4023 ifirst=max(ipixstart,(ixp2-1)*nxif1+ixpmin1)
4024 ilast=min(ipixend,(ixp2-1)*nxif1+ixpmax1)
4025 if (ifirst>ilast) cycle
4026 do ipix=ifirst,ilast
4027 ixp1=1+mod(ipix-1,nxif1)
4029 profile_batch(1)=profile_batch(1)+one
4031 cache(iigrid)%box_max,hit,t_enter,t_exit)
4033 profile_batch(2)=profile_batch(2)+one
4036 cache(iigrid)%opacity,cache(iigrid)%sourcev,&
4037 ipix,ray_origin,cache(iigrid)%xface1,&
4038 cache(iigrid)%xface2,cache(iigrid)%xface3,&
4040 segments,nseg,capacity)
4041 profile_batch(3)=profile_batch(3)+dble(nseg-nsegbefore)
4047 call mpi_allreduce(nseg,maxnsegbatch,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
4048 if (maxnsegbatch>maxsegbatchtarget .and. npixbatch>1)
then
4049 npixbatch=max(1,npixbatch/2)
4050 ipixend=ipixstart+npixbatch-1
4051 if (
allocated(segments))
deallocate(segments)
4054 batchaccepted=.true.
4058 profile_local=profile_local+profile_batch
4060 write(*,
'(a,3(i0,1x))')
' cart_dda thick adaptive batch: ',&
4061 ipixstart,ipixend,maxnsegbatch
4064 if (.not.
allocated(segments))
then
4066 allocate(segments(nsegvars,capacity))
4071 ownersegcounts(owner)=ownersegcounts(owner)+1
4073 sendcounts=nsegvars*ownersegcounts
4076 senddispls(ipe)=senddispls(ipe-1)+sendcounts(ipe-1)
4079 allocate(segments_send(nsegvars,max(1,nseg)))
4083 isegdest=senddispls(owner)/nsegvars+owneroffsets(owner)+1
4084 segments_send(:,isegdest)=segments(:,is)
4085 owneroffsets(owner)=owneroffsets(owner)+1
4088 call mpi_alltoall(sendcounts,1,mpi_integer,recvcounts,1,mpi_integer,
icomm,
ierrmpi)
4091 recvdispls(ipe)=recvdispls(ipe-1)+recvcounts(ipe-1)
4093 totalcount=sum(recvcounts)
4094 totalseg=totalcount/nsegvars
4095 profile_local(4)=profile_local(4)+dble(totalcount)
4096 allocate(segments_recv(nsegvars,max(1,totalseg)))
4099 maxownersegcountlocal=maxval(ownersegcounts)
4100 call mpi_allreduce(maxownersegcountlocal,maxownersegcount,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
4101 do segoffset=0,maxownersegcount-1,maxsegcommtarget
4103 roundsenddispls=senddispls
4105 if (ownersegcounts(ipe)>segoffset)
then
4106 roundsendcounts(ipe)=nsegvars*min(maxsegcommtarget,ownersegcounts(ipe)-segoffset)
4107 roundsenddispls(ipe)=senddispls(ipe)+nsegvars*segoffset
4111 call mpi_alltoall(roundsendcounts,1,mpi_integer,roundrecvcounts,1,mpi_integer,
icomm,
ierrmpi)
4112 roundrecvdispls(0)=0
4114 roundrecvdispls(ipe)=roundrecvdispls(ipe-1)+roundrecvcounts(ipe-1)
4116 totalroundcount=sum(roundrecvcounts)
4117 totalroundseg=totalroundcount/nsegvars
4118 allocate(segments_recv_round(nsegvars,max(1,totalroundseg)))
4120 call mpi_alltoallv(segments_send,roundsendcounts,roundsenddispls,mpi_double_precision,&
4121 segments_recv_round,roundrecvcounts,roundrecvdispls,&
4124 if (totalroundseg>0)
then
4125 segments_recv(:,recvfill+1:recvfill+totalroundseg)=segments_recv_round(:,1:totalroundseg)
4126 recvfill=recvfill+totalroundseg
4128 deallocate(segments_recv_round)
4131 if (recvfill/=totalseg)
call mpistop(
"cart_dda thick segmented receive mismatch")
4133 if (totalseg>0)
then
4134 allocate(idx(totalseg))
4135 bucketcounts(1:npixbatch)=0
4138 ipix=nint(segments_recv(1,is))
4140 ilocal=ipix-ipixstart+1
4141 bucketcounts(ilocal)=bucketcounts(ilocal)+1
4147 do ilocal=1,npixbatch
4148 bucketoffsets(ilocal+1)=bucketoffsets(ilocal)+bucketcounts(ilocal)
4150 bucketfill(1:npixbatch)=bucketoffsets(1:npixbatch)
4153 ipix=nint(segments_recv(1,is))
4155 ilocal=ipix-ipixstart+1
4156 idx(bucketfill(ilocal))=is
4157 bucketfill(ilocal)=bucketfill(ilocal)+1
4162 do ipix=ipixstart,ipixend
4164 ilocal=ipix-ipixstart+1
4165 nidx=bucketcounts(ilocal)
4167 profile_local(5)=profile_local(5)+dble(nidx)*dble(nidx)
4169 ixglobal=1+mod(ipix-1,nxif1)
4170 iyglobal=1+(ipix-1)/nxif1
4171 do iseg=bucketoffsets(ilocal),bucketoffsets(ilocal+1)-1
4173 euvthin(ixglobal,iyglobal)=euvthin(ixglobal,iyglobal)+segments_recv(3,is)
4175 euv(ixglobal,iyglobal)=euv(ixglobal,iyglobal)+atten*segments_recv(3,is)
4176 dpl(ixglobal,iyglobal)=dpl(ixglobal,iyglobal)+atten*segments_recv(5,is)
4177 tau(ixglobal,iyglobal)=tau(ixglobal,iyglobal)+max(zero,segments_recv(4,is))
4184 deallocate(segments_send,segments_recv)
4185 if (
allocated(segments))
deallocate(segments)
4189 do iigrid=1,igridstail
4190 if (
allocated(cache(iigrid)%source))
deallocate(cache(iigrid)%source)
4191 if (
allocated(cache(iigrid)%opacity))
deallocate(cache(iigrid)%opacity)
4192 if (
allocated(cache(iigrid)%sourcev))
deallocate(cache(iigrid)%sourcev)
4193 if (
allocated(cache(iigrid)%xface1))
deallocate(cache(iigrid)%xface1)
4194 if (
allocated(cache(iigrid)%xface2))
deallocate(cache(iigrid)%xface2)
4195 if (
allocated(cache(iigrid)%xface3))
deallocate(cache(iigrid)%xface3)
4198 deallocate(sendcounts,recvcounts,senddispls,recvdispls,roundsendcounts,roundrecvcounts,&
4199 roundsenddispls,roundrecvdispls,ownersegcounts,owneroffsets,bucketcounts,&
4200 bucketoffsets,bucketfill)
4201 allocate(image_reduce(nxif1,nxif2))
4202 call mpi_allreduce(euv,image_reduce,nxif1*nxif2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4204 call mpi_allreduce(dpl,image_reduce,nxif1*nxif2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4206 call mpi_allreduce(tau,image_reduce,nxif1*nxif2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4208 call mpi_allreduce(euvthin,image_reduce,nxif1*nxif2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4209 euvthin=image_reduce
4210 deallocate(image_reduce)
4211 call mpi_allreduce(profile_local,profile_global,5,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4213 write(*,
'(a,5(es12.5,1x))') &
4214 ' cart_dda thick profile: ',profile_global
4221 integer,
intent(in) :: nXIF1,nXIF2
4223 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
4224 double precision,
intent(out) :: Tau(nXIF1,nXIF2),EUVthin(nXIF1,nXIF2)
4226 integer :: ixO^L,ixO^D,ixI^L,ix^D
4227 integer :: iigrid,igrid,levelg,rft,direction_LOS,nLOS,numSeg,nLayerVars,nLayerSeg
4228 integer :: ixP1,ixP2,ixL,iSub1,iSub2,relL
4229 integer :: nLosBatch,nBatch,iBatch,ixLstart,ixLend,ixLgridStart,ixLgridEnd
4230 integer :: ixPmin1,ixPmin2
4231 double precision :: ds_cm,jds,kds,jvds,atten,layerBytes,targetBytes
4232 double precision,
allocatable :: rho(:^D&)
4233 double precision,
allocatable :: layer_ds(:,:,:,:),layer_all(:,:,:,:)
4248 if (nxif1>huge(numseg)/max(1,nxif2) .or. nxif1*nxif2>huge(numseg)/nlayervars)
then
4249 call mpistop(
"thick EUV layer buffer is too large for one MPI reduction")
4251 nlayerseg=nxif1*nxif2*nlayervars
4252 targetbytes=256.d0*1024.d0*1024.d0
4253 layerbytes=dble(nlayerseg)*8.d0*2.d0
4254 nlosbatch=max(1,min(16,int(targetbytes/max(one,layerbytes))))
4255 if (nlayerseg>huge(numseg)/nlosbatch)
then
4256 call mpistop(
"thick EUV batched layer buffer is too large for one MPI reduction")
4259 allocate(cache(igridstail))
4260 do iigrid=1,igridstail; igrid=igrids(iigrid);
4261 ^d&ixomin^d=ixmlo^d\
4262 ^d&ixomax^d=ixmhi^d\
4263 ^d&iximin^d=
ixglo^d\
4264 ^d&iximax^d=
ixghi^d\
4266 cache(iigrid)%igrid=igrid
4267 levelg=ps(igrid)%level
4269 cache(iigrid)%level=levelg
4270 cache(iigrid)%rft=rft
4272 select case(direction_los)
4274 cache(iigrid)%los_min=(
node(pig1_,igrid)-1)*rft*block_nx1+1
4275 cache(iigrid)%los_max=
node(pig1_,igrid)*rft*block_nx1
4277 cache(iigrid)%los_min=(
node(pig2_,igrid)-1)*rft*block_nx2+1
4278 cache(iigrid)%los_max=
node(pig2_,igrid)*rft*block_nx2
4280 cache(iigrid)%los_min=(
node(pig3_,igrid)-1)*rft*block_nx3+1
4281 cache(iigrid)%los_max=
node(pig3_,igrid)*rft*block_nx3
4284 allocate(cache(iigrid)%source(ixi^s),cache(iigrid)%opacity(ixi^s),&
4285 cache(iigrid)%sourcev(ixi^s),rho(ixi^s))
4286 cache(iigrid)%source=zero
4287 cache(iigrid)%opacity=zero
4288 cache(iigrid)%sourcev=zero
4291 cache(iigrid)%source,cache(iigrid)%opacity)
4293 call get_euv(
wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,cache(iigrid)%source)
4296 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
4297 cache(iigrid)%sourcev(ixo^s)=cache(iigrid)%source(ixo^s)*&
4298 (-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s))
4304 allocate(layer_ds(nxif1,nxif2,nlayervars,nlosbatch),layer_all(nxif1,nxif2,nlayervars,nlosbatch))
4310 do ixlstart=1,nlos,nlosbatch
4311 ixlend=min(nlos,ixlstart+nlosbatch-1)
4312 nbatch=ixlend-ixlstart+1
4313 layer_ds(:,:,:,1:nbatch)=zero
4315 do iigrid=1,igridstail
4316 ixlgridstart=max(ixlstart,cache(iigrid)%los_min)
4317 ixlgridend=min(ixlend,cache(iigrid)%los_max)
4318 if (ixlgridstart>ixlgridend) cycle
4319 igrid=cache(iigrid)%igrid
4320 rft=cache(iigrid)%rft
4321 ^d&ixomin^d=ixmlo^d\
4322 ^d&ixomax^d=ixmhi^d\
4323 ^d&iximin^d=
ixglo^d\
4324 ^d&iximax^d=
ixghi^d\
4326 do ixl=ixlgridstart,ixlgridend
4327 ibatch=ixl-ixlstart+1
4328 rell=ixl-cache(iigrid)%los_min
4330 select case(direction_los)
4332 ix1=ixomin1+rell/rft
4333 do ix2=ixomin2,ixomax2
4334 ixpmin1=(
node(pig2_,igrid)-1)*rft*block_nx2+(ix2-ixomin2)*rft+1
4335 do ix3=ixomin3,ixomax3
4336 ixpmin2=(
node(pig3_,igrid)-1)*rft*block_nx3+(ix3-ixomin3)*rft+1
4339 jds=cache(iigrid)%source(ix^d)*ds_cm
4340 kds=cache(iigrid)%opacity(ix^d)*ds_cm
4341 jvds=cache(iigrid)%sourcev(ix^d)*ds_cm
4346 layer_ds(ixp1,ixp2,1,ibatch)=layer_ds(ixp1,ixp2,1,ibatch)+jds
4347 layer_ds(ixp1,ixp2,2,ibatch)=layer_ds(ixp1,ixp2,2,ibatch)+kds
4348 layer_ds(ixp1,ixp2,3,ibatch)=layer_ds(ixp1,ixp2,3,ibatch)+jvds
4354 ix2=ixomin2+rell/rft
4355 do ix3=ixomin3,ixomax3
4356 ixpmin1=(
node(pig3_,igrid)-1)*rft*block_nx3+(ix3-ixomin3)*rft+1
4357 do ix1=ixomin1,ixomax1
4358 ixpmin2=(
node(pig1_,igrid)-1)*rft*block_nx1+(ix1-ixomin1)*rft+1
4361 jds=cache(iigrid)%source(ix^d)*ds_cm
4362 kds=cache(iigrid)%opacity(ix^d)*ds_cm
4363 jvds=cache(iigrid)%sourcev(ix^d)*ds_cm
4368 layer_ds(ixp1,ixp2,1,ibatch)=layer_ds(ixp1,ixp2,1,ibatch)+jds
4369 layer_ds(ixp1,ixp2,2,ibatch)=layer_ds(ixp1,ixp2,2,ibatch)+kds
4370 layer_ds(ixp1,ixp2,3,ibatch)=layer_ds(ixp1,ixp2,3,ibatch)+jvds
4376 ix3=ixomin3+rell/rft
4377 do ix1=ixomin1,ixomax1
4378 ixpmin1=(
node(pig1_,igrid)-1)*rft*block_nx1+(ix1-ixomin1)*rft+1
4379 do ix2=ixomin2,ixomax2
4380 ixpmin2=(
node(pig2_,igrid)-1)*rft*block_nx2+(ix2-ixomin2)*rft+1
4383 jds=cache(iigrid)%source(ix^d)*ds_cm
4384 kds=cache(iigrid)%opacity(ix^d)*ds_cm
4385 jvds=cache(iigrid)%sourcev(ix^d)*ds_cm
4390 layer_ds(ixp1,ixp2,1,ibatch)=layer_ds(ixp1,ixp2,1,ibatch)+jds
4391 layer_ds(ixp1,ixp2,2,ibatch)=layer_ds(ixp1,ixp2,2,ibatch)+kds
4392 layer_ds(ixp1,ixp2,3,ibatch)=layer_ds(ixp1,ixp2,3,ibatch)+jvds
4401 numseg=nlayerseg*nbatch
4402 call mpi_allreduce(layer_ds,layer_all,numseg,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4408 euvthin(ixp1,ixp2)=euvthin(ixp1,ixp2)+layer_all(ixp1,ixp2,1,ibatch)
4410 euv(ixp1,ixp2)=euv(ixp1,ixp2)+atten*layer_all(ixp1,ixp2,1,ibatch)
4411 dpl(ixp1,ixp2)=dpl(ixp1,ixp2)+atten*layer_all(ixp1,ixp2,3,ibatch)
4412 tau(ixp1,ixp2)=tau(ixp1,ixp2)+layer_all(ixp1,ixp2,2,ibatch)
4419 do iigrid=1,igridstail
4420 if (
allocated(cache(iigrid)%source))
deallocate(cache(iigrid)%source)
4421 if (
allocated(cache(iigrid)%opacity))
deallocate(cache(iigrid)%opacity)
4422 if (
allocated(cache(iigrid)%sourcev))
deallocate(cache(iigrid)%sourcev)
4424 deallocate(cache,layer_ds,layer_all)
4433 integer,
intent(in) :: ixI^L, ixO^L
4434 double precision,
intent(in) :: x(ixI^S,1:ndim),dx(ixI^S,1:ndim)
4435 double precision,
allocatable,
intent(out) :: rface(:),thetaface(:),phiface(:)
4436 integer :: ix1,ix2,ix3
4438 allocate(rface(ixomin1:ixomax1+1),thetaface(ixomin2:ixomax2+1),phiface(ixomin3:ixomax3+1))
4439 do ix1=ixomin1,ixomax1
4440 rface(ix1)=x(ix1,ixomin2,ixomin3,1)-half*dx(ix1,ixomin2,ixomin3,1)
4442 rface(ixomax1+1)=x(ixomax1,ixomin2,ixomin3,1)+half*dx(ixomax1,ixomin2,ixomin3,1)
4443 do ix2=ixomin2,ixomax2
4444 thetaface(ix2)=x(ixomin1,ix2,ixomin3,2)-half*dx(ixomin1,ix2,ixomin3,2)
4446 thetaface(ixomax2+1)=x(ixomin1,ixomax2,ixomin3,2)+half*dx(ixomin1,ixomax2,ixomin3,2)
4447 do ix3=ixomin3,ixomax3
4448 phiface(ix3)=x(ixomin1,ixomin2,ix3,3)-half*dx(ixomin1,ixomin2,ix3,3)
4450 phiface(ixomax3+1)=x(ixomin1,ixomin2,ixomax3,3)+half*dx(ixomin1,ixomin2,ixomax3,3)
4454 double precision,
allocatable,
intent(inout) :: tvals(:)
4455 integer,
intent(inout) :: nt,capacity
4456 double precision,
intent(in) :: t
4458 double precision,
allocatable :: tmp(:)
4460 if (t /= t .or. abs(t)>1.d90)
return
4461 if (.not.
allocated(tvals))
then
4463 allocate(tvals(capacity))
4464 else if (nt>=capacity)
then
4465 allocate(tmp(capacity))
4468 allocate(tvals(2*capacity))
4469 tvals(1:capacity)=tmp
4478 double precision,
intent(inout) :: tvals(:)
4479 integer,
intent(inout) :: nt
4480 double precision,
intent(in) :: t
4482 if (t /= t .or. abs(t)>1.d90)
return
4483 if (nt>=
size(tvals))
return
4489 double precision,
intent(inout) :: tvals(:)
4490 integer,
intent(inout) :: nt
4493 double precision :: key,epsT
4499 do while (j>=1 .and. tvals(j)>key)
4505 epst=max(1.d-12,1.d-10*max(one,abs(tvals(nt)-tvals(1))))
4508 if (abs(tvals(i)-tvals(nout))>epst)
then
4510 tvals(nout)=tvals(i)
4517 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),rface
4518 double precision,
allocatable,
intent(inout) :: tvals(:)
4519 integer,
intent(inout) :: nt,capacity
4521 double precision :: aa,bb,cc,disc,root
4524 bb=2.d0*sum(ray_origin*ray_dir)
4525 cc=sum(ray_origin**2)-rface**2
4526 disc=bb**2-4.d0*aa*cc
4527 if (disc<zero)
return
4528 root=sqrt(max(zero,disc))
4529 call sph_add_t(tvals,nt,capacity,(-bb-root)/(2.d0*aa))
4530 call sph_add_t(tvals,nt,capacity,(-bb+root)/(2.d0*aa))
4534 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),thetaface
4535 double precision,
allocatable,
intent(inout) :: tvals(:)
4536 integer,
intent(inout) :: nt,capacity
4538 double precision :: cth,aa,bb,cc,disc,root
4544 if (abs(cth)<1.d-12)
then
4545 if (abs(ray_dir(3))>1.d-14) &
4546 call sph_add_t(tvals,nt,capacity,-ray_origin(3)/ray_dir(3))
4549 aa=ray_dir(3)**2-cth**2*sum(ray_dir**2)
4550 bb=2.d0*(ray_origin(3)*ray_dir(3)-cth**2*sum(ray_origin*ray_dir))
4551 cc=ray_origin(3)**2-cth**2*sum(ray_origin**2)
4552 if (abs(aa)<1.d-14)
then
4553 if (abs(bb)>1.d-14)
call sph_add_t(tvals,nt,capacity,-cc/bb)
4556 disc=bb**2-4.d0*aa*cc
4557 if (disc<zero)
return
4558 root=sqrt(max(zero,disc))
4559 call sph_add_t(tvals,nt,capacity,(-bb-root)/(2.d0*aa))
4560 call sph_add_t(tvals,nt,capacity,(-bb+root)/(2.d0*aa))
4564 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),phiface
4565 double precision,
allocatable,
intent(inout) :: tvals(:)
4566 integer,
intent(inout) :: nt,capacity
4568 double precision :: normal(1:3),denom,numer
4570 normal(1)=-sin(phiface)
4571 normal(2)=cos(phiface)
4573 denom=sum(normal*ray_dir)
4574 if (abs(denom)<1.d-14)
return
4575 numer=sum(normal*ray_origin)
4576 call sph_add_t(tvals,nt,capacity,-numer/denom)
4580 double precision,
intent(in) :: pos(1:3)
4581 double precision,
intent(out) :: sph(1:3)
4583 sph(1)=sqrt(sum(pos**2))
4584 if (sph(1)>zero)
then
4585 sph(2)=acos(max(-one,min(one,pos(3)/sph(1))))
4589 sph(3)=atan2(pos(2),pos(1))
4593 integer,
intent(in) :: imin,imax
4594 double precision,
intent(in) ::
value,faces(imin:imax+1)
4596 integer :: ilo,ihi,imid
4599 if (
value<faces(imin)-1.d-12 .or.
value>faces(imax+1)+1.d-12)
return
4600 if (
value<=faces(imin))
then
4604 if (
value>=faces(imax+1))
then
4611 do while (ihi-ilo>1)
4613 if (
value>=faces(imid))
then
4619 idx=min(imax,max(imin,ilo))
4623 integer,
intent(in) :: imin,imax
4624 double precision,
intent(in) ::
value,faces(imin:imax+1)
4626 integer :: ilo,ihi,imid
4629 if (
value>faces(imin)+1.d-12 .or.
value<faces(imax+1)-1.d-12)
return
4630 if (
value>=faces(imin))
then
4634 if (
value<=faces(imax+1))
then
4641 do while (ihi-ilo>1)
4643 if (
value<=faces(imid))
then
4649 idx=min(imax,max(imin,ilo))
4653 double precision,
intent(in) :: pos(1:3)
4654 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4655 phiface(ixOmin3:ixOmax3+1)
4656 integer,
intent(in) :: ixO^L
4657 integer,
intent(out) :: ix1,ix2,ix3
4658 logical,
intent(out) :: inside
4660 double precision :: sph(1:3),phi
4664 if (phi<phiface(ixomin3)-1.d-12) phi=phi+2.d0*dpi
4665 if (phi>phiface(ixomax3+1)+1.d-12) phi=phi-2.d0*dpi
4669 inside=ix1>0 .and. ix2>0 .and. ix3>0
4674 double precision,
intent(in) :: pos(1:3),ximg1,ximg2
4676 double precision :: dotp,rc,rthick,rloc
4679 rc=sqrt(ximg1**2+ximg2**2)
4680 rloc=sqrt(sum(pos**2))
4683 if (dotp>=
zero)
then
4684 if (rc<=rthick) visible=.false.
4686 if (rloc<=rthick) visible=.false.
4691 ixPmin1,ixPmax1,ixPmin2,ixPmax2,has_pixels)
4692 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4693 phiface(ixOmin3:ixOmax3+1)
4694 integer,
intent(in) :: ixO^L,nXI1,nXI2
4695 double precision,
intent(in) :: xI1(nXI1),xI2(nXI2),dxI
4696 integer,
intent(out) :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
4697 logical,
intent(out) :: has_pixels
4699 integer,
parameter :: nsample=5
4701 double precision :: sph(1:3),xcent(1:2)
4702 double precision :: xmin1,xmax1,xmin2,xmax2
4703 double precision :: wr,wt,wp,pad
4711 wr=dble(ir)/dble(nsample-1)
4712 sph(1)=(one-wr)*rface(ixomin1)+wr*rface(ixomax1+1)
4714 wt=dble(it)/dble(nsample-1)
4715 sph(2)=(one-wt)*thetaface(ixomin2)+wt*thetaface(ixomax2+1)
4717 wp=dble(ip)/dble(nsample-1)
4718 if (ir/=0 .and. ir/=nsample-1 .and. it/=0 .and. it/=nsample-1 .and. &
4719 ip/=0 .and. ip/=nsample-1) cycle
4720 sph(3)=(one-wp)*phiface(ixomin3)+wp*phiface(ixomax3+1)
4722 xmin1=min(xmin1,xcent(1))
4723 xmax1=max(xmax1,xcent(1))
4724 xmin2=min(xmin2,xcent(2))
4725 xmax2=max(xmax2,xcent(2))
4734 ixpmin1=max(1,floor((xmin1-(xi1(1)-half*dxi))/dxi)+1)
4735 ixpmax1=min(nxi1,ceiling((xmax1-(xi1(1)-half*dxi))/dxi))
4736 ixpmin2=max(1,floor((xmin2-(xi2(1)-half*dxi))/dxi)+1)
4737 ixpmax2=min(nxi2,ceiling((xmax2-(xi2(1)-half*dxi))/dxi))
4738 has_pixels=ixpmin1<=ixpmax1 .and. ixpmin2<=ixpmax2
4742 rface,thetaface,phiface,EUVp)
4743 integer,
intent(in) :: ixI^L,ixO^L
4744 double precision,
intent(in) :: source(ixI^S),ray_origin(1:3),ximg1,ximg2
4745 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4746 phiface(ixOmin3:ixOmax3+1)
4747 double precision,
intent(inout) :: EUVp
4749 integer :: nt,capacity,i,ix^D
4750 double precision,
allocatable :: tvals(:)
4751 double precision :: posMid(1:3),ds_cm,tMid,t0,t1
4756 do ix1=ixomin1,ixomax1+1
4759 do ix2=ixomin2,ixomax2+1
4762 do ix3=ixomin3,ixomax3+1
4766 if (
allocated(tvals))
deallocate(tvals)
4776 posmid=ray_origin+tmid*
vec_los
4778 call sph_locate_cell(posmid,rface,thetaface,phiface,ixo^l,ix1,ix2,ix3,inside)
4779 if (.not. inside) cycle
4780 ds_cm=(t1-t0)*unit_length
4781 if (si_unit) ds_cm=ds_cm*1.d2
4782 euvp=euvp+source(ix^d)*ds_cm
4788 pixel_id,ray_origin,ximg1,ximg2,&
4789 rface,thetaface,phiface,rface2,&
4790 theta_cos,phi_sin,phi_cos,&
4791 segments,nseg,capacity)
4794 integer,
intent(in) :: ixI^L,ixO^L,pixel_id
4795 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S)
4796 double precision,
intent(in) :: ray_origin(1:3),ximg1,ximg2
4797 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4798 phiface(ixOmin3:ixOmax3+1)
4799 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4800 phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
4801 double precision,
allocatable,
intent(inout) :: segments(:,:)
4802 integer,
intent(inout) :: nseg,capacity
4804 integer :: nt,i,ix^D
4805 double precision :: tvals(2*(ixOmax1-ixOmin1+2)+2*(ixOmax2-ixOmin2+2)+&
4806 (ixOmax3-ixOmin3+2))
4807 double precision :: posMid(1:3),ds_cm,tMid,t0,t1,jds,kds
4808 double precision :: dir2,origin2,odotd,aa,bb,cc,disc,root,cth2
4809 double precision :: denom,numer,r2,mu,phi,dotp,rthick2,rc2
4814 origin2=sum(ray_origin**2)
4816 rthick2=(r_opt_thick*
const_rsun/unit_length)**2
4817 rc2=ximg1**2+ximg2**2
4819 do ix1=ixomin1,ixomax1+1
4822 cc=origin2-rface2(ix1)
4823 disc=bb**2-4.d0*aa*cc
4824 if (disc>=
zero)
then
4825 root=sqrt(max(
zero,disc))
4830 do ix2=ixomin2,ixomax2+1
4831 if (abs(theta_cos(ix2))<1.d-12)
then
4836 cth2=theta_cos(ix2)**2
4838 bb=2.d0*(ray_origin(3)*
vec_los(3)-cth2*odotd)
4839 cc=ray_origin(3)**2-cth2*origin2
4840 if (abs(aa)<1.d-14)
then
4843 disc=bb**2-4.d0*aa*cc
4844 if (disc>=
zero)
then
4845 root=sqrt(max(
zero,disc))
4851 do ix3=ixomin3,ixomax3+1
4853 if (abs(denom)>=1.d-14)
then
4854 numer=-phi_sin(ix3)*ray_origin(1)+phi_cos(ix3)*ray_origin(2)
4866 posmid=ray_origin+tmid*
vec_los
4869 if (dotp>=
zero)
then
4870 if (rc2<=rthick2) cycle
4872 if (r2<=rthick2) cycle
4877 mu=posmid(3)/sqrt(r2)
4882 phi=atan2(posmid(2),posmid(1))
4883 if (phi<phiface(ixomin3)-1.d-12) phi=phi+2.d0*
dpi
4884 if (phi>phiface(ixomax3+1)+1.d-12) phi=phi-2.d0*
dpi
4886 inside=ix1>0 .and. ix2>0 .and. ix3>0
4887 if (.not. inside) cycle
4888 ds_cm=(t1-t0)*unit_length
4889 if (si_unit) ds_cm=ds_cm*1.d2
4890 jds=max(
zero,source(ix^d))*ds_cm
4891 kds=max(
zero,opacity(ix^d))*ds_cm
4897 double precision,
intent(in) :: pos(1:3)
4898 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4899 phiface(ixOmin3:ixOmax3+1)
4900 integer,
intent(in) :: ixO^L
4901 integer,
intent(out) :: ix1,ix2,ix3
4902 logical,
intent(out) :: inside
4904 double precision :: r2,mu,phi
4914 phi=atan2(pos(2),pos(1))
4915 if (phi<phiface(ixomin3)-1.d-12) phi=phi+2.d0*dpi
4916 if (phi>phiface(ixomax3+1)+1.d-12) phi=phi-2.d0*dpi
4918 inside=ix1>0 .and. ix2>0 .and. ix3>0
4922 double precision,
intent(in) :: t,tNow,tExit,epsRay
4923 double precision,
intent(inout) :: tNext
4924 logical,
intent(inout) :: found
4926 if (t>tnow+epsray .and. t<=texit+epsray .and. t<tnext)
then
4933 double precision,
intent(in) :: t,theta_face_cos,ray_origin(1:3),tNow,tExit,epsRay
4934 double precision,
intent(inout) :: tNext
4935 logical,
intent(inout) :: found
4937 double precision :: pos(1:3),r2
4939 if (t<=tnow+epsray .or. t>texit+epsray .or. t>=tnext)
return
4942 if (r2<=zero)
return
4943 if (theta_face_cos>1.d-12 .and. pos(3)<-1.d-10)
return
4944 if (theta_face_cos<-1.d-12 .and. pos(3)>1.d-10)
return
4945 if (abs(pos(3)**2-theta_face_cos**2*r2)>1.d-6*max(one,r2))
return
4951 double precision,
intent(in) :: t,phi_face_sin,phi_face_cos,ray_origin(1:3),tNow,tExit,epsRay
4952 double precision,
intent(inout) :: tNext
4953 logical,
intent(inout) :: found
4955 double precision :: pos(1:3),radialDot
4957 if (t<=tnow+epsray .or. t>texit+epsray .or. t>=tnext)
return
4959 radialdot=phi_face_cos*pos(1)+phi_face_sin*pos(2)
4960 if (radialdot<-1.d-10)
return
4966 ix1,ix2,ix3,tNow,tExit,epsRay,tNext,found)
4967 double precision,
intent(in) :: ray_origin(1:3)
4968 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4969 phiface(ixOmin3:ixOmax3+1)
4970 double precision,
intent(in) :: phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
4971 integer,
intent(in) :: ixO^L,ix1,ix2,ix3
4972 double precision,
intent(in) :: tNow,tExit,epsRay
4973 double precision,
intent(out) :: tNext
4974 logical,
intent(out) :: found
4977 double precision :: dir2,origin2,odotd,aa,bb,cc,disc,root,cth2,denom,numer
4982 origin2=sum(ray_origin**2)
4988 cc=origin2-rface2(iface)
4989 disc=bb**2-4.d0*aa*cc
4990 if (disc>=zero)
then
4991 root=sqrt(max(zero,disc))
4998 if (abs(theta_cos(iface))<1.d-12)
then
5000 -ray_origin(3)/
vec_los(3),theta_cos(iface),ray_origin,tnow,texit,epsray,&
5004 cth2=theta_cos(iface)**2
5006 bb=2.d0*(ray_origin(3)*
vec_los(3)-cth2*odotd)
5007 cc=ray_origin(3)**2-cth2*origin2
5008 if (abs(aa)<1.d-14)
then
5010 ray_origin,tnow,texit,epsray,tnext,found)
5012 disc=bb**2-4.d0*aa*cc
5013 if (disc>=zero)
then
5014 root=sqrt(max(zero,disc))
5016 ray_origin,tnow,texit,epsray,tnext,found)
5018 ray_origin,tnow,texit,epsray,tnext,found)
5025 if (abs(denom)>=1.d-14)
then
5026 numer=-phi_sin(iface)*ray_origin(1)+phi_cos(iface)*ray_origin(2)
5028 tnow,texit,epsray,tnext,found)
5034 ray_origin,ximg1,ximg2,rface2,theta_cos,phiface,&
5036 t_enter,t_exit,segments,nseg,capacity,ok)
5039 integer,
intent(in) :: ixI^L,ixO^L,pixel_id
5040 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S)
5041 double precision,
intent(in) :: ray_origin(1:3),ximg1,ximg2
5042 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
5043 phiface(ixOmin3:ixOmax3+1)
5044 double precision,
intent(in) :: phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
5045 double precision,
intent(in) :: t_enter,t_exit
5046 double precision,
allocatable,
intent(inout) :: segments(:,:)
5047 integer,
intent(inout) :: nseg,capacity
5048 logical,
intent(out) :: ok
5050 integer :: ix^D,nstep,maxSteps
5051 double precision :: tNow,tNext,tEnd,tMid,epsRay,ds_cm,jds,kds
5052 double precision :: pos(1:3),r2,dotp,rthick2,rc2
5053 logical :: inside,found
5056 if (t_exit<=t_enter)
return
5057 epsray=max(1.d-12,1.d-10*max(
one,abs(t_exit-t_enter)))
5058 pos=ray_origin+(t_enter+epsray)*
vec_los
5060 if (.not. inside)
then
5065 rthick2=(r_opt_thick*
const_rsun/unit_length)**2
5066 rc2=ximg1**2+ximg2**2
5069 maxsteps=8*((ixomax1-ixomin1+1)+(ixomax2-ixomin2+1)+(ixomax3-ixomin3+1)+3)
5071 do while (tnow<t_exit-epsray)
5073 ix1,ix2,ix3,tnow,t_exit,epsray,tnext,found)
5074 if (.not. found)
then
5078 tend=min(tnext,t_exit)
5080 tmid=
half*(tnow+tend)
5084 if (.not. ((dotp>=
zero .and. rc2<=rthick2) .or. &
5085 (dotp<
zero .and. r2<=rthick2)))
then
5086 ds_cm=(tend-tnow)*unit_length
5087 if (si_unit) ds_cm=ds_cm*1.d2
5088 jds=max(
zero,source(ix^d))*ds_cm
5089 kds=max(
zero,opacity(ix^d))*ds_cm
5094 if (tnow>=t_exit-epsray)
exit
5095 pos=ray_origin+(tnow+epsray)*
vec_los
5097 if (.not. inside)
then
5102 if (nstep>maxsteps)
then
5110 pixel_id,ray_origin,ximg1,ximg2,&
5111 rface,thetaface,phiface,rface2,&
5112 theta_cos,phi_sin,phi_cos,&
5113 segments,nseg,capacity,fallback)
5114 integer,
intent(in) :: ixI^L,ixO^L,pixel_id
5115 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S)
5116 double precision,
intent(in) :: ray_origin(1:3),ximg1,ximg2
5117 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
5118 phiface(ixOmin3:ixOmax3+1)
5119 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
5120 phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
5121 double precision,
allocatable,
intent(inout) :: segments(:,:)
5122 integer,
intent(inout) :: nseg,capacity
5123 logical,
intent(out) :: fallback
5125 integer :: nt,i,nsegStart,ix^D
5126 double precision :: tvals(12),posMid(1:3),t0,t1,tMid
5127 double precision :: dir2,origin2,odotd,aa,bb,cc,disc,root,cth2,denom,numer
5128 logical :: inside,ok
5134 origin2=sum(ray_origin**2)
5137 do ix1=ixomin1,ixomax1+1,ixomax1-ixomin1+1
5140 cc=origin2-rface2(ix1)
5141 disc=bb**2-4.d0*aa*cc
5142 if (disc>=zero)
then
5143 root=sqrt(max(zero,disc))
5148 do ix2=ixomin2,ixomax2+1,ixomax2-ixomin2+1
5149 if (abs(theta_cos(ix2))<1.d-12)
then
5154 cth2=theta_cos(ix2)**2
5156 bb=2.d0*(ray_origin(3)*
vec_los(3)-cth2*odotd)
5157 cc=ray_origin(3)**2-cth2*origin2
5158 if (abs(aa)<1.d-14)
then
5161 disc=bb**2-4.d0*aa*cc
5162 if (disc>=zero)
then
5163 root=sqrt(max(zero,disc))
5169 do ix3=ixomin3,ixomax3+1,ixomax3-ixomin3+1
5171 if (abs(denom)>=1.d-14)
then
5172 numer=-phi_sin(ix3)*ray_origin(1)+phi_cos(ix3)*ray_origin(2)
5185 posmid=ray_origin+tmid*
vec_los
5187 if (.not. inside) cycle
5189 ray_origin,ximg1,ximg2,rface2,theta_cos,phiface,phi_sin,phi_cos,&
5190 t0,t1,segments,nseg,capacity,ok)
5195 pixel_id,ray_origin,ximg1,ximg2,rface,thetaface,phiface,rface2,&
5196 theta_cos,phi_sin,phi_cos,segments,nseg,capacity)
5206 integer,
intent(in) :: numXI1,numXI2
5207 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2),dxI
5209 double precision,
intent(inout) :: EM(numXI1,numXI2)
5211 integer :: ixO^L,ixI^L,ix^D
5212 integer :: iigrid,igrid,ixP1,ixP2,ixPmin1,ixPmax1,ixPmin2,ixPmax2
5213 integer :: iseg,nseg,capacity,sphDdaFallbackLocal,sphDdaFallbackGlobal
5214 double precision :: ray_origin(1:3),profile_local(3),profile_global(3)
5215 double precision,
allocatable :: source(:^D&)
5216 double precision,
allocatable :: rface(:),thetaface(:),phiface(:)
5217 double precision,
allocatable :: rface2(:),theta_cos(:),phi_sin(:),phi_cos(:)
5218 double precision,
allocatable :: segments(:,:)
5219 logical :: has_pixels,ddaFallback
5222 sphddafallbacklocal=0
5223 do iigrid=1,igridstail; igrid=igrids(iigrid);
5224 ^d&ixomin^d=ixmlo^d\
5225 ^d&ixomax^d=ixmhi^d\
5226 ^d&iximin^d=
ixglo^d\
5227 ^d&iximax^d=
ixghi^d\
5229 allocate(source(ixi^s))
5234 allocate(rface2(ixomin1:ixomax1+1),theta_cos(ixomin2:ixomax2+1),&
5235 phi_sin(ixomin3:ixomax3+1),phi_cos(ixomin3:ixomax3+1))
5237 theta_cos=cos(thetaface)
5238 phi_sin=sin(phiface)
5239 phi_cos=cos(phiface)
5242 ixpmin1,ixpmax1,ixpmin2,ixpmax2,has_pixels)
5243 if (has_pixels)
then
5245 do ixp1=ixpmin1,ixpmax1
5246 do ixp2=ixpmin2,ixpmax2
5248 profile_local(1)=profile_local(1)+one
5252 1,ray_origin,xi1(ixp1),xi2(ixp2),rface,thetaface,phiface,&
5253 rface2,theta_cos,phi_sin,phi_cos,segments,nseg,capacity,ddafallback)
5254 if (ddafallback) sphddafallbacklocal=sphddafallbacklocal+1
5256 em(ixp1,ixp2)=em(ixp1,ixp2)+segments(3,iseg)
5260 rface,thetaface,phiface,em(ixp1,ixp2))
5264 profile_local(2)=profile_local(2)+dble((ixpmax1-ixpmin1+1)*(ixpmax2-ixpmin2+1))
5266 if (
allocated(segments))
deallocate(segments)
5267 profile_local(3)=profile_local(3)+one
5268 if (
allocated(rface2))
deallocate(rface2)
5269 if (
allocated(theta_cos))
deallocate(theta_cos)
5270 if (
allocated(phi_sin))
deallocate(phi_sin)
5271 if (
allocated(phi_cos))
deallocate(phi_cos)
5272 deallocate(source,rface,thetaface,phiface)
5274 call mpi_allreduce(profile_local,profile_global,3,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5275 call mpi_allreduce(sphddafallbacklocal,sphddafallbackglobal,1,mpi_integer,mpi_sum,
icomm,
ierrmpi)
5278 write(*,
'(a,3(es12.5,1x))')
' sph_dda thin profile rays pixels blocks: ',profile_global
5279 write(*,
'(a,i0)')
' sph_dda thin fallback rays: ',sphddafallbackglobal
5281 write(*,
'(a,3(es12.5,1x))')
' sph_intersection thin profile rays pixels blocks: ',profile_global
5289 integer,
intent(in) :: numXI1,numXI2
5290 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2),dxI
5292 double precision,
intent(out) :: EUV(numXI1,numXI2),Tau(numXI1,numXI2),EUVthin(numXI1,numXI2)
5294 integer,
parameter :: nSegVars=5
5295 integer :: ixO^L,ixI^L,ix^D
5296 integer :: iigrid,igrid,ixP1,ixP2,ipix,ipixStart,ipixEnd,nPixBatch,pixel_id
5297 integer :: nseg,capacity,totalCount,totalSeg,ipe,is,iseg,nidx,owner,isegDest,nsegBefore
5298 integer :: ixGlobal,iyGlobal,ixPmin1,ixPmax1,ixPmin2,ixPmax2,iFirst,iLast,iLocal
5299 integer :: nPixBatchTarget
5300 integer :: maxSegBatchTarget,maxSegCommTarget,maxNsegBatch,nPixTotal
5301 integer :: maxOwnerSegCount,maxOwnerSegCountLocal,segOffset,recvFill,totalRoundCount,totalRoundSeg
5302 integer :: sphDdaFallbackLocal,sphDdaFallbackGlobal
5303 integer,
allocatable :: sendCounts(:),recvCounts(:),sendDispls(:),recvDispls(:)
5304 integer,
allocatable :: roundSendCounts(:),roundRecvCounts(:)
5305 integer,
allocatable :: roundSendDispls(:),roundRecvDispls(:)
5306 integer,
allocatable :: ownerSegCounts(:),ownerOffsets(:),idx(:)
5307 integer,
allocatable :: bucketCounts(:),bucketOffsets(:),bucketFill(:)
5308 double precision :: ray_origin(1:3),atten
5309 double precision :: profile_local(5),profile_global(5),profile_batch(5)
5310 double precision :: phys_max_local(2),phys_max_global(2)
5311 double precision :: phys_sum_local(2),phys_sum_global(2),phys_sum_batch(2)
5312 double precision,
allocatable :: segments(:,:),segments_send(:,:),segments_recv(:,:)
5313 double precision,
allocatable :: segments_recv_round(:,:)
5314 double precision,
allocatable :: image_reduce(:,:)
5315 logical :: has_pixels,batchAccepted,batchReduced,ddaFallback
5324 sphddafallbacklocal=0
5325 allocate(sendcounts(0:
npe-1),recvcounts(0:
npe-1),senddispls(0:
npe-1),recvdispls(0:
npe-1))
5326 allocate(roundsendcounts(0:
npe-1),roundrecvcounts(0:
npe-1))
5327 allocate(roundsenddispls(0:
npe-1),roundrecvdispls(0:
npe-1))
5328 allocate(ownersegcounts(0:
npe-1),owneroffsets(0:
npe-1))
5329 allocate(cache(igridstail))
5331 allocate(bucketcounts(npixbatchtarget),bucketoffsets(npixbatchtarget+1),&
5332 bucketfill(npixbatchtarget))
5334 do iigrid=1,igridstail; igrid=igrids(iigrid);
5335 ^d&ixomin^d=ixmlo^d\
5336 ^d&ixomax^d=ixmhi^d\
5337 ^d&iximin^d=
ixglo^d\
5338 ^d&iximax^d=
ixghi^d\
5340 cache(iigrid)%igrid=igrid
5341 allocate(cache(iigrid)%source(ixi^s),cache(iigrid)%opacity(ixi^s))
5342 cache(iigrid)%source=zero
5343 cache(iigrid)%opacity=zero
5344 call get_euv(
wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,cache(iigrid)%source)
5347 phys_max_local(1)=max(phys_max_local(1),maxval(cache(iigrid)%source(ixo^s)))
5348 phys_max_local(2)=max(phys_max_local(2),maxval(cache(iigrid)%opacity(ixo^s)))
5350 cache(iigrid)%rface,cache(iigrid)%thetaface,&
5351 cache(iigrid)%phiface)
5352 allocate(cache(iigrid)%rface2(ixomin1:ixomax1+1),&
5353 cache(iigrid)%theta_cos(ixomin2:ixomax2+1),&
5354 cache(iigrid)%phi_sin(ixomin3:ixomax3+1),&
5355 cache(iigrid)%phi_cos(ixomin3:ixomax3+1))
5356 cache(iigrid)%rface2=cache(iigrid)%rface**2
5357 cache(iigrid)%theta_cos=cos(cache(iigrid)%thetaface)
5358 cache(iigrid)%phi_sin=sin(cache(iigrid)%phiface)
5359 cache(iigrid)%phi_cos=cos(cache(iigrid)%phiface)
5361 cache(iigrid)%phiface,ixo^l,numxi1,numxi2,xi1,xi2,dxi,&
5362 cache(iigrid)%ixPmin1,cache(iigrid)%ixPmax1,&
5363 cache(iigrid)%ixPmin2,cache(iigrid)%ixPmax2,cache(iigrid)%has_pixels)
5366 npixtotal=numxi1*numxi2
5368 do while (ipixstart<=npixtotal)
5369 ipixend=min(numxi1*numxi2,ipixstart+npixbatchtarget-1)
5370 npixbatch=ipixend-ipixstart+1
5371 batchaccepted=.false.
5372 batchreduced=.false.
5374 do while (.not. batchaccepted)
5380 do iigrid=1,igridstail; igrid=igrids(iigrid);
5381 ^d&ixomin^d=ixmlo^d\
5382 ^d&ixomax^d=ixmhi^d\
5383 ^d&iximin^d=
ixglo^d\
5384 ^d&iximax^d=
ixghi^d\
5386 ixpmin1=cache(iigrid)%ixPmin1
5387 ixpmax1=cache(iigrid)%ixPmax1
5388 ixpmin2=cache(iigrid)%ixPmin2
5389 ixpmax2=cache(iigrid)%ixPmax2
5390 has_pixels=cache(iigrid)%has_pixels
5391 if (.not. has_pixels) cycle
5393 do ixp2=ixpmin2,ixpmax2
5394 ifirst=max(ipixstart,(ixp2-1)*numxi1+ixpmin1)
5395 ilast=min(ipixend,(ixp2-1)*numxi1+ixpmax1)
5396 if (ifirst>ilast) cycle
5397 do ipix=ifirst,ilast
5398 ixp1=1+mod(ipix-1,numxi1)
5401 profile_batch(1)=profile_batch(1)+one
5405 cache(iigrid)%opacity,pixel_id,ray_origin,xi1(ixp1),xi2(ixp2),&
5406 cache(iigrid)%rface,cache(iigrid)%thetaface,cache(iigrid)%phiface,&
5407 cache(iigrid)%rface2,cache(iigrid)%theta_cos,&
5408 cache(iigrid)%phi_sin,cache(iigrid)%phi_cos,&
5409 segments,nseg,capacity,ddafallback)
5410 if (ddafallback) sphddafallbacklocal=sphddafallbacklocal+1
5413 cache(iigrid)%opacity,pixel_id,ray_origin,xi1(ixp1),xi2(ixp2),&
5414 cache(iigrid)%rface,cache(iigrid)%thetaface,cache(iigrid)%phiface,&
5415 cache(iigrid)%rface2,cache(iigrid)%theta_cos,&
5416 cache(iigrid)%phi_sin,cache(iigrid)%phi_cos,&
5417 segments,nseg,capacity)
5419 if (nseg>nsegbefore) profile_batch(2)=profile_batch(2)+one
5420 profile_batch(3)=profile_batch(3)+dble(nseg-nsegbefore)
5421 do iseg=nsegbefore+1,nseg
5422 phys_sum_batch(1)=phys_sum_batch(1)+segments(3,iseg)
5423 phys_sum_batch(2)=phys_sum_batch(2)+segments(4,iseg)
5429 call mpi_allreduce(nseg,maxnsegbatch,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
5430 if (maxnsegbatch>maxsegbatchtarget .and. npixbatch>1)
then
5431 npixbatch=max(1,npixbatch/2)
5432 ipixend=ipixstart+npixbatch-1
5433 if (
allocated(segments))
deallocate(segments)
5436 batchaccepted=.true.
5440 profile_local=profile_local+profile_batch
5441 phys_sum_local=phys_sum_local+phys_sum_batch
5444 write(*,
'(a,3(i0,1x))')
' sph_dda thick adaptive batch: ',&
5445 ipixstart,ipixend,maxnsegbatch
5447 write(*,
'(a,3(i0,1x))')
' sph_intersection thick adaptive batch: ',&
5448 ipixstart,ipixend,maxnsegbatch
5452 if (.not.
allocated(segments))
then
5454 allocate(segments(nsegvars,capacity))
5459 ownersegcounts(owner)=ownersegcounts(owner)+1
5461 sendcounts=nsegvars*ownersegcounts
5464 senddispls(ipe)=senddispls(ipe-1)+sendcounts(ipe-1)
5467 allocate(segments_send(nsegvars,max(1,nseg)))
5471 isegdest=senddispls(owner)/nsegvars+owneroffsets(owner)+1
5472 segments_send(:,isegdest)=segments(:,is)
5473 owneroffsets(owner)=owneroffsets(owner)+1
5476 call mpi_alltoall(sendcounts,1,mpi_integer,recvcounts,1,mpi_integer,
icomm,
ierrmpi)
5479 recvdispls(ipe)=recvdispls(ipe-1)+recvcounts(ipe-1)
5481 totalcount=sum(recvcounts)
5482 totalseg=totalcount/nsegvars
5483 profile_local(4)=profile_local(4)+dble(totalcount)
5484 allocate(segments_recv(nsegvars,max(1,totalseg)))
5487 maxownersegcountlocal=maxval(ownersegcounts)
5488 call mpi_allreduce(maxownersegcountlocal,maxownersegcount,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
5489 do segoffset=0,maxownersegcount-1,maxsegcommtarget
5491 roundsenddispls=senddispls
5493 if (ownersegcounts(ipe)>segoffset)
then
5494 roundsendcounts(ipe)=nsegvars*min(maxsegcommtarget,ownersegcounts(ipe)-segoffset)
5495 roundsenddispls(ipe)=senddispls(ipe)+nsegvars*segoffset
5499 call mpi_alltoall(roundsendcounts,1,mpi_integer,roundrecvcounts,1,mpi_integer,
icomm,
ierrmpi)
5500 roundrecvdispls(0)=0
5502 roundrecvdispls(ipe)=roundrecvdispls(ipe-1)+roundrecvcounts(ipe-1)
5504 totalroundcount=sum(roundrecvcounts)
5505 totalroundseg=totalroundcount/nsegvars
5506 allocate(segments_recv_round(nsegvars,max(1,totalroundseg)))
5508 call mpi_alltoallv(segments_send,roundsendcounts,roundsenddispls,mpi_double_precision,&
5509 segments_recv_round,roundrecvcounts,roundrecvdispls,&
5512 if (totalroundseg>0)
then
5513 segments_recv(:,recvfill+1:recvfill+totalroundseg)=segments_recv_round(:,1:totalroundseg)
5514 recvfill=recvfill+totalroundseg
5516 deallocate(segments_recv_round)
5519 if (recvfill/=totalseg)
call mpistop(
"ray-segment receive mismatch")
5521 if (totalseg>0)
then
5522 allocate(idx(totalseg))
5523 bucketcounts(1:npixbatch)=0
5526 ipix=nint(segments_recv(1,is))
5528 ilocal=ipix-ipixstart+1
5529 bucketcounts(ilocal)=bucketcounts(ilocal)+1
5535 do ilocal=1,npixbatch
5536 bucketoffsets(ilocal+1)=bucketoffsets(ilocal)+bucketcounts(ilocal)
5538 bucketfill(1:npixbatch)=bucketoffsets(1:npixbatch)
5541 ipix=nint(segments_recv(1,is))
5543 ilocal=ipix-ipixstart+1
5544 idx(bucketfill(ilocal))=is
5545 bucketfill(ilocal)=bucketfill(ilocal)+1
5550 do ipix=ipixstart,ipixend
5552 ilocal=ipix-ipixstart+1
5553 nidx=bucketcounts(ilocal)
5555 profile_local(5)=profile_local(5)+dble(nidx)*dble(nidx)
5557 idx(bucketoffsets(ilocal):bucketoffsets(ilocal+1)-1),nidx)
5558 ixglobal=1+mod(ipix-1,numxi1)
5559 iyglobal=1+(ipix-1)/numxi1
5560 do iseg=bucketoffsets(ilocal),bucketoffsets(ilocal+1)-1
5562 euvthin(ixglobal,iyglobal)=euvthin(ixglobal,iyglobal)+segments_recv(3,is)
5564 euv(ixglobal,iyglobal)=euv(ixglobal,iyglobal)+atten*segments_recv(3,is)
5565 tau(ixglobal,iyglobal)=tau(ixglobal,iyglobal)+max(zero,segments_recv(4,is))
5572 deallocate(segments_send,segments_recv)
5573 if (
allocated(segments))
deallocate(segments)
5577 do iigrid=1,igridstail
5578 if (
allocated(cache(iigrid)%source))
deallocate(cache(iigrid)%source)
5579 if (
allocated(cache(iigrid)%opacity))
deallocate(cache(iigrid)%opacity)
5580 if (
allocated(cache(iigrid)%rface))
deallocate(cache(iigrid)%rface)
5581 if (
allocated(cache(iigrid)%thetaface))
deallocate(cache(iigrid)%thetaface)
5582 if (
allocated(cache(iigrid)%phiface))
deallocate(cache(iigrid)%phiface)
5583 if (
allocated(cache(iigrid)%rface2))
deallocate(cache(iigrid)%rface2)
5584 if (
allocated(cache(iigrid)%theta_cos))
deallocate(cache(iigrid)%theta_cos)
5585 if (
allocated(cache(iigrid)%phi_sin))
deallocate(cache(iigrid)%phi_sin)
5586 if (
allocated(cache(iigrid)%phi_cos))
deallocate(cache(iigrid)%phi_cos)
5589 deallocate(sendcounts,recvcounts,senddispls,recvdispls,roundsendcounts,roundrecvcounts,&
5590 roundsenddispls,roundrecvdispls,ownersegcounts,owneroffsets,bucketcounts,&
5591 bucketoffsets,bucketfill)
5592 allocate(image_reduce(numxi1,numxi2))
5593 call mpi_allreduce(euv,image_reduce,numxi1*numxi2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5595 call mpi_allreduce(tau,image_reduce,numxi1*numxi2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5597 call mpi_allreduce(euvthin,image_reduce,numxi1*numxi2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5598 euvthin=image_reduce
5599 deallocate(image_reduce)
5600 call mpi_allreduce(profile_local,profile_global,5,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5601 call mpi_allreduce(phys_max_local,phys_max_global,2,mpi_double_precision,mpi_max,
icomm,
ierrmpi)
5602 call mpi_allreduce(phys_sum_local,phys_sum_global,2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5603 call mpi_allreduce(sphddafallbacklocal,sphddafallbackglobal,1,mpi_integer,mpi_sum,
icomm,
ierrmpi)
5606 write(*,
'(a,5(es12.5,1x))')
' sph_dda thick profile: ',profile_global
5607 write(*,
'(a,i0)')
' sph_dda thick fallback rays: ',sphddafallbackglobal
5608 write(*,
'(a,4(es12.5,1x))')
' sph_dda thick physics maxj maxk sumjds sumkds: ',&
5609 phys_max_global(1),phys_max_global(2),phys_sum_global(1),phys_sum_global(2)
5611 write(*,
'(a,5(es12.5,1x))')
' sph_intersection thick profile: ',profile_global
5612 write(*,
'(a,4(es12.5,1x))')
' sph_intersection thick physics maxj maxk sumjds sumkds: ',&
5613 phys_max_global(1),phys_max_global(2),phys_sum_global(1),phys_sum_global(2)
5623 double precision,
intent(out) :: xImin1,xImax1,xImin2,xImax2
5625 integer,
parameter :: nsample=5
5626 integer :: iigrid,igrid,ir,it,ip
5627 integer :: ixI^L,ixO^L
5628 double precision,
allocatable :: rface(:),thetaface(:),phiface(:)
5629 double precision :: local_min1,local_max1,local_min2,local_max2
5630 double precision :: sph(1:3),xcent(1:2),wr,wt,wp
5632 local_min1=huge(one)
5633 local_max1=-huge(one)
5634 local_min2=huge(one)
5635 local_max2=-huge(one)
5637 do iigrid=1,igridstail
5638 igrid=igrids(iigrid)
5639 ^d&ixomin^d=ixmlo^d\
5640 ^d&ixomax^d=ixmhi^d\
5641 ^d&iximin^d=ixglo^d\
5642 ^d&iximax^d=ixghi^d\
5645 rface,thetaface,phiface)
5647 wr=dble(ir)/dble(nsample-1)
5648 sph(1)=(one-wr)*rface(ixomin1)+wr*rface(ixomax1+1)
5650 wt=dble(it)/dble(nsample-1)
5651 sph(2)=(one-wt)*thetaface(ixomin2)+wt*thetaface(ixomax2+1)
5653 wp=dble(ip)/dble(nsample-1)
5654 if (ir/=0 .and. ir/=nsample-1 .and. it/=0 .and. it/=nsample-1 .and. &
5655 ip/=0 .and. ip/=nsample-1) cycle
5656 sph(3)=(one-wp)*phiface(ixomin3)+wp*phiface(ixomax3+1)
5658 local_min1=min(local_min1,xcent(1))
5659 local_max1=max(local_max1,xcent(1))
5660 local_min2=min(local_min2,xcent(2))
5661 local_max2=max(local_max2,xcent(2))
5665 deallocate(rface,thetaface,phiface)
5668 call mpi_allreduce(local_min1,ximin1,1,mpi_double_precision,mpi_min,icomm,ierrmpi)
5669 call mpi_allreduce(local_max1,ximax1,1,mpi_double_precision,mpi_max,icomm,ierrmpi)
5670 call mpi_allreduce(local_min2,ximin2,1,mpi_double_precision,mpi_min,icomm,ierrmpi)
5671 call mpi_allreduce(local_max2,ximax2,1,mpi_double_precision,mpi_max,icomm,ierrmpi)
5672 if (ximin1>0.5d0*huge(one) .or. ximax1<-0.5d0*huge(one) .or. &
5673 ximin2>0.5d0*huge(one) .or. ximax2<-0.5d0*huge(one))
then
5674 call mpistop(
"sph_intersection could not determine image bounds")
5679 double precision,
intent(out) :: dxI
5681 select case(trim(dat_resolution_mode))
5687 call mpistop(
"unknown dat_resolution_mode")
5692 double precision,
intent(out) :: dxI
5694 double precision :: refine_factor,dr,dtheta,dphi,rmin,sin_theta_min
5696 refine_factor=dble(2**(refine_max_level-1))
5698 dxi=min(abs(xprobmax1-xprobmin1)/(dble(domain_nx1)*refine_factor),&
5699 abs(xprobmax2-xprobmin2)/(dble(domain_nx2)*refine_factor),&
5700 abs(xprobmax3-xprobmin3)/(dble(domain_nx3)*refine_factor))
5702 rmin=max(smalldouble,min(xprobmin1,xprobmax1))
5703 sin_theta_min=max(smalldouble,min(abs(sin(xprobmin2)),&
5704 abs(sin(xprobmax2))))
5705 dr=abs(xprobmax1-xprobmin1)/(dble(domain_nx1)*refine_factor)
5706 dtheta=abs(xprobmax2-xprobmin2)/(dble(domain_nx2)*refine_factor)
5707 dphi=abs(xprobmax3-xprobmin3)/(dble(domain_nx3)*refine_factor)
5708 dxi=min(dr,rmin*dtheta,rmin*sin_theta_min*dphi)
5710 call mpistop(
"nominal dat resolution needs Cartesian or spherical coordinates")
5713 if (dxi<=zero .or. dxi>half*huge(one))
then
5714 call mpistop(
"could not determine nominal dat-resolution image spacing")
5719 double precision,
intent(out) :: dxI
5721 integer :: iigrid,igrid,ixI^L,ixO^L,ix^D
5722 double precision :: local_min,global_min,dr,ds_theta,ds_phi,rval,theta
5725 do iigrid=1,igridstail
5726 igrid=igrids(iigrid)
5727 ^d&ixomin^d=ixmlo^d\
5728 ^d&ixomax^d=ixmhi^d\
5729 ^d&iximin^d=ixglo^d\
5730 ^d&iximax^d=ixghi^d\
5732 do ix1=ixomin1,ixomax1
5733 do ix2=ixomin2,ixomax2
5734 do ix3=ixomin3,ixomax3
5736 local_min=min(local_min,ps(igrid)%dx(ix^d,1),&
5737 ps(igrid)%dx(ix^d,2),ps(igrid)%dx(ix^d,3))
5739 rval=max(smalldouble,ps(igrid)%x(ix^d,1))
5740 theta=ps(igrid)%x(ix^d,2)
5741 dr=ps(igrid)%dx(ix^d,1)
5742 ds_theta=rval*ps(igrid)%dx(ix^d,2)
5743 ds_phi=rval*max(smalldouble,sin(theta))*ps(igrid)%dx(ix^d,3)
5744 local_min=min(local_min,dr,ds_theta,ds_phi)
5746 call mpistop(
"minimum dat resolution needs Cartesian or spherical coordinates")
5753 call mpi_allreduce(local_min,global_min,1,mpi_double_precision,mpi_min,&
5755 if (global_min<=zero .or. global_min>half*huge(one))
then
5756 call mpistop(
"could not determine minimum dat-resolution image spacing")
5767 integer,
intent(in) :: qunit
5769 character(20),
intent(in) :: datatype
5771 integer :: ix^D,numXI1,numXI2,numWI
5772 double precision :: xImin1,xImax1,xImin2,xImax2,xIcent1,xIcent2,dxI
5773 double precision,
allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:)
5774 double precision,
allocatable :: wI(:,:,:),wIs(:,:,:),EM(:,:),Dpl(:,:),Tau(:,:),EMthin(:,:),WLB(:,:,:)
5775 double precision :: vec_temp1(1:3),vec_temp2(1:3)
5776 double precision :: vec_z(1:3),vec_cor(1:3),xI_cor(1:2)
5777 double precision :: res,LOS_psi,r_max,r_loc
5780 character (30) :: ion
5781 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
5782 double precision :: arcsec,RHESSI_rsl,LASCO_rsl,pixel,R_occult,smallflux
5783 integer :: iigrid,igrid,i,j,numSI,iw
5784 logical :: emit,ray_image_global,has_thick_output
5798 ximin1=-abs(xprobmax1)
5799 ximin2=-abs(xprobmax1)
5800 ximax1=abs(xprobmax1)
5801 ximax2=abs(xprobmax1)
5806 if (ix1==1) vec_cor(1)=xprobmin1
5807 if (ix1==2) vec_cor(1)=xprobmax1
5809 if (ix2==1) vec_cor(2)=xprobmin2
5810 if (ix2==2) vec_cor(2)=xprobmax2
5812 if (ix3==1) vec_cor(3)=xprobmin3
5813 if (ix3==2) vec_cor(3)=xprobmax3
5816 r_loc=r_loc+(vec_cor(2)-
x_origin(2))**2
5817 r_loc=r_loc+(vec_cor(3)-
x_origin(3))**2
5819 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
5822 r_max=max(r_max,r_loc)
5826 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
5832 ximin1=min(ximin1,xi_cor(1))
5833 ximax1=max(ximax1,xi_cor(1))
5834 ximin2=min(ximin2,xi_cor(2))
5835 ximax2=max(ximax2,xi_cor(2))
5848 xicent1=(ximin1+ximax1)/2.d0
5849 xicent2=(ximin2+ximax2)/2.d0
5857 if (datatype==
'image_euv')
then
5861 else if (datatype==
'image_sxr')
then
5863 dxi=rhessi_rsl*arcsec
5865 else if (datatype==
'image_whitelight')
then
5876 call mpistop(
'Whitelight synthesis: instrument is not supported!')
5878 dxi=lasco_rsl*arcsec
5883 numxi1=8*ceiling((ximax1-xicent1)/dxi/8.d0)
5884 ximin1=xicent1-numxi1*dxi
5885 ximax1=xicent1+numxi1*dxi
5887 numxi2=8*ceiling((ximax2-xicent2)/dxi/8.d0)
5888 ximin2=xicent2-numxi2*dxi
5889 ximax2=xicent2+numxi2*dxi
5891 allocate(xi1(numxi1),xi2(numxi2),dxi1(numxi1),dxi2(numxi2))
5893 xi1(ix1)=ximin1+dxi*(ix1-
half)
5897 xi2(ix2)=ximin2+dxi*(ix2-
half)
5902 if (datatype==
'image_euv' .or. datatype==
'image_sxr')
then
5903 has_thick_output=datatype==
'image_euv' .and. trim(
radiation_transfer)==
'thick' .and. &
5906 if (datatype==
'image_euv')
then
5911 allocate(wi(numxi1,numxi2,numwi),wis(numxi1,numxi2,numwi),em(numxi1,numxi2))
5915 ray_image_global=.false.
5916 if (has_thick_output)
then
5917 allocate(tau(numxi1,numxi2),emthin(numxi1,numxi2))
5921 if (
slab .and. datatype==
'image_euv' .and. &
5923 ray_image_global=.true.
5924 allocate(dpl(numxi1,numxi2))
5933 do iigrid=1,igridstail; igrid=igrids(iigrid);
5936 else if (trim(
ray_method_active) ==
'spherical' .and. datatype ==
'image_euv')
then
5938 ray_image_global=.true.
5944 do iigrid=1,igridstail; igrid=igrids(iigrid);
5948 if (ray_image_global)
then
5949 if (has_thick_output)
then
5951 has_thick_output,tau=tau,euvthin=emthin,&
5952 cap_absorption=.true.)
5959 if (em(ix1,ix2)>smallflux) wis(ix1,ix2,1)=em(ix1,ix2)
5963 if (.not. ray_image_global)
then
5964 numsi=numxi1*numxi2*numwi
5965 call mpi_allreduce(wis,wi,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5973 call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
5974 if (
allocated(tau))
deallocate(tau)
5975 if (
allocated(emthin))
deallocate(emthin)
5976 deallocate(wi,wis,em)
5977 else if (datatype==
'image_whitelight')
then
5979 allocate(wi(numxi1,numxi2,numwi),wis(numxi1,numxi2,numwi),wlb(numxi1,numxi2,numwi))
5984 do iigrid=1,igridstail; igrid=igrids(iigrid);
5990 if (wlb(ix1,ix2,1)>smallflux)
then
5991 wis(ix1,ix2,1)=wlb(ix1,ix2,1)
5992 wis(ix1,ix2,2)=wlb(ix1,ix2,2)
5996 numsi=numxi1*numxi2*numwi
5997 call mpi_allreduce(wis,wi,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
6004 call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
6005 deallocate(wi,wis,wlb)
6008 deallocate(xi1,xi2,dxi1,dxi2)
6013 integer,
intent(in) :: igrid,numXI1,numXI2
6014 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
6015 double precision,
intent(in) :: dxI
6017 character(20),
intent(in) :: datatype
6018 double precision,
intent(inout) :: EM(numXI1,numXI2)
6020 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
6021 double precision :: xb^L,xd^D
6022 double precision,
allocatable :: flux(:^D&),opacity(:^D&)
6023 double precision :: res
6024 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
6025 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC
6026 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2)
6029 double precision :: logTe
6030 character (30) :: ion
6031 double precision :: lineCent
6032 double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
6033 double precision :: arcsec,pixel,RHESSI_rsl,area_1AU
6034 double precision :: aa,bb
6036 ^d&ixomin^d=ixmlo^d\
6037 ^d&ixomax^d=ixmhi^d\
6038 ^d&iximin^d=ixglo^d\
6039 ^d&iximax^d=ixghi^d\
6040 ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
6041 ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
6044 arcsec=7.25d5/unit_length
6046 arcsec=7.25d7/unit_length
6049 allocate(flux(ixi^s),opacity(ixi^s))
6050 if (datatype==
'image_euv')
then
6051 if (trim(emission_model)==
'pseudo_current')
then
6053 else if (trim(emission_model)==
'radio_ff')
then
6057 call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
6058 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
6060 call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
6061 pixel=spacersl*arcsec
6062 sigma0=sigma_psf*pixel
6063 else if (datatype==
'image_sxr')
then
6065 call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr)
6066 rhessi_rsl=2.3d0/instrument_resolution_factor
6068 pixel=rhessi_rsl*arcsec
6069 sigma0=sigma_psf*pixel
6074 {
do ix^d=ixomin^d,ixomax^d\}
6076 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi1(^d))/(dxi/2.d0)));
6077 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi2(^d))/(dxi/2.d0)));
6078 ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
6079 if (datatype==
'image_euv')
then
6081 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length*1.d2/dxi/dxi
6083 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length/dxi/dxi
6085 else if (datatype==
'image_sxr')
then
6087 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length**3/area_1au
6089 if (fluxsubc>smalldouble)
then
6091 {
do isubc^d=1,nsubc^d\}
6092 ^d&xsubc(^d)=ps(igrid)%x(ix^dd,^d)-half*ps(igrid)%dx(ix^dd,^d)+(isubc^d-half)*dxsubc^d;
6096 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
6097 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
6098 ixpmin1=max(1,ixp1-3)
6099 ixpmax1=min(ixp1+3,numxi1)
6100 ixpmin2=max(1,ixp2-3)
6101 ixpmax2=min(ixp2+3,numxi2)
6102 do ixp1=ixpmin1,ixpmax1
6103 do ixp2=ixpmin2,ixpmax2
6104 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6105 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6106 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6107 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6108 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
6109 em(ixp1,ixp2)=em(ixp1,ixp2)+fluxsubc*factor
6116 deallocate(flux,opacity)
6120 integer,
intent(in) :: igrid,numXI1,numXI2
6121 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
6122 double precision,
intent(in) :: dxI
6124 character(20),
intent(in) :: datatype
6125 double precision,
intent(inout) :: EM(numXI1,numXI2)
6127 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
6128 double precision,
allocatable :: flux(:^D&),Ne(:^D&),opacity(:^D&)
6129 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
6130 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC,RsubC
6131 double precision :: TBsubC,PBsubC
6132 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2),xSubC_car(1:3)
6133 double precision :: R_thick,dotp,dvolume,R_occult,Rc
6134 double precision :: dxl(1:3),x_sph(1:3),dx_sph(1:3)
6135 double precision :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
6136 logical :: sun_back_side,emit
6139 double precision :: logTe
6140 character (30) :: ion
6141 double precision :: lineCent
6142 double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
6143 double precision :: RHESSI_rsl,area_1AU,arcsec,pixel
6145 ^d&ixomin^d=ixmlo^d;
6146 ^d&ixomax^d=ixmhi^d;
6147 ^d&iximin^d=ixglo^d;
6148 ^d&iximax^d=ixghi^d;
6151 arcsec=7.25d5/unit_length
6153 arcsec=7.25d7/unit_length
6156 allocate(flux(ixi^s),opacity(ixi^s))
6157 if (datatype==
'image_euv')
then
6158 if (trim(emission_model)==
'pseudo_current')
then
6160 else if (trim(emission_model)==
'radio_ff')
then
6164 call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
6165 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
6167 call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
6168 pixel=spacersl*arcsec
6169 sigma0=sigma_psf*pixel
6170 else if (datatype==
'image_sxr')
then
6172 call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr)
6173 rhessi_rsl=2.3d0/instrument_resolution_factor
6175 pixel=rhessi_rsl*arcsec
6176 sigma0=sigma_psf*pixel
6181 r_thick=r_opt_thick*const_rsun/unit_length
6182 {
do ix^d=ixomin^d,ixomax^d\}
6183 x_sph(1:3)=ps(igrid)%x(ix^d,1:3)
6184 dx_sph(1:3)=ps(igrid)%dx(ix^d,1:3)
6186 dxl(2)=x_sph(1)*dx_sph(2)
6187 dxl(3)=x_sph(1)*dsin(x_sph(2))*dx_sph(3)
6192 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
6194 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
6196 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
6198 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
6200 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
6202 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
6207 xsubc(1)=x_sph(1)-half*dx_sph(1)+(isubc1-half)*dx_sph(1)/nsubc1
6209 dxsubc1=dx_sph(1)/nsubc1
6212 xsubc(2)=x_sph(2)-half*dx_sph(2)+(isubc2-half)*dx_sph(2)/nsubc2
6213 dxsubc2=xsubc(1)*dx_sph(2)/nsubc2
6214 dxsubc3=xsubc(1)*dsin(xsubc(2))*dx_sph(3)/nsubc3
6215 dvolume=dxsubc1*dxsubc2*dxsubc3
6216 if (datatype==
'image_euv')
then
6218 fluxsubc=flux(ix^d)*dvolume*unit_length*1.d2/dxi/dxi
6220 fluxsubc=flux(ix^d)*dvolume*unit_length/dxi/dxi
6222 else if (datatype==
'image_sxr')
then
6224 fluxsubc=flux(ix^d)*dvolume*unit_length**3/area_1au
6227 if (fluxsubc>smalldouble)
then
6230 xsubc(3)=x_sph(3)-half*dx_sph(3)+(isubc3-half)*dx_sph(3)/nsubc3
6232 rc=dsqrt(xcent(1)**2+xcent(2)**2)
6237 sun_back_side=.true.
6238 if (dotp<0.d0) sun_back_side=.false.
6240 if (sun_back_side)
then
6242 if (rc>r_thick) emit=.true.
6245 if (xsubc(1)<=r_thick) emit=.false.
6251 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
6252 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
6253 ixpmin1=max(1,ixp1-3)
6254 ixpmax1=min(ixp1+3,numxi1)
6255 ixpmin2=max(1,ixp2-3)
6256 ixpmax2=min(ixp2+3,numxi2)
6257 do ixp1=ixpmin1,ixpmax1
6258 do ixp2=ixpmin2,ixpmax2
6259 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6260 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6261 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6262 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6263 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
6264 em(ixp1,ixp2)=em(ixp1,ixp2)+fluxsubc*factor
6274 deallocate(flux,opacity)
6281 integer,
intent(in) :: igrid,numXI1,numXI2,numWI
6282 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
6283 double precision,
intent(in) :: dxI
6285 character(20),
intent(in) :: datatype
6286 double precision,
intent(inout) :: WLB(numXI1,numXI2,numWI)
6288 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
6289 double precision,
allocatable :: flux(:^D&),Ne(:^D&),nH_dummy(:^D&)
6290 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
6291 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC,RsubC
6292 double precision :: sigma_PSF,sigma0,arcsec,pixel,LASCO_rsl
6293 double precision :: A,B,C,D,Rc,Ne0,TBsubC,PBsubC,factor
6294 double precision :: R_thick,dotp,dvolume,R_occult
6295 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2),xSubC_car(1:3)
6296 double precision :: dxl(1:3),x_sph(1:3),dx_sph(1:3)
6297 double precision :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
6300 ^d&ixomin^d=ixmlo^d;
6301 ^d&ixomax^d=ixmhi^d;
6302 ^d&iximin^d=ixglo^d;
6303 ^d&iximax^d=ixghi^d;
6306 arcsec=7.25d5/unit_length
6308 arcsec=7.25d7/unit_length
6312 allocate(nh_dummy(ixi^s))
6313 if (whitelight_instrument==
'LASCO/C1')
then
6314 lasco_rsl=5.6d0/instrument_resolution_factor
6316 else if (whitelight_instrument==
'LASCO/C2')
then
6317 lasco_rsl=11.4d0/instrument_resolution_factor
6319 else if (whitelight_instrument==
'LASCO/C3')
then
6320 lasco_rsl=56.d0/instrument_resolution_factor
6323 if (r_occultor>1.d0) r_occult=r_occultor
6324 r_occult=r_occult*const_rsun/unit_length
6325 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,ne)
6327 call eos%get_ne_nH(ixi^l, ixo^l, ps(igrid)%w, ps(igrid)%x, ne, nh_dummy)
6329 pixel=lasco_rsl*arcsec
6330 sigma0=sigma_psf*pixel
6333 r_thick=r_opt_thick*const_rsun/unit_length
6334 {
do ix^d=ixomin^d,ixomax^d\}
6335 x_sph(1:3)=ps(igrid)%x(ix^d,1:3)
6336 dx_sph(1:3)=ps(igrid)%dx(ix^d,1:3)
6338 dxl(2)=x_sph(1)*dx_sph(2)
6339 dxl(3)=x_sph(1)*dsin(x_sph(2))*dx_sph(3)
6340 ne0=ne(ix^d)*unit_numberdensity
6345 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
6347 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
6349 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
6351 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
6353 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
6355 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
6360 xsubc(1)=x_sph(1)-half*dx_sph(1)+(isubc1-half)*dx_sph(1)/nsubc1
6362 dxsubc1=dx_sph(1)/nsubc1
6366 xsubc(2)=x_sph(2)-half*dx_sph(2)+(isubc2-half)*dx_sph(2)/nsubc2
6367 dxsubc2=xsubc(1)*dx_sph(2)/nsubc2
6368 dxsubc3=xsubc(1)*dsin(xsubc(2))*dx_sph(3)/nsubc3
6369 dvolume=dxsubc1*dxsubc2*dxsubc3
6372 xsubc(3)=x_sph(3)-half*dx_sph(3)+(isubc3-half)*dx_sph(3)/nsubc3
6374 rc=dsqrt(xcent(1)**2+xcent(2)**2)
6377 if (rc>r_occult)
then
6381 tbsubc=tbsubc*dvolume*unit_length/dxi/dxi
6382 pbsubc=pbsubc*dvolume*unit_length/dxi/dxi
6383 if (tbsubc<1.d-20) emit=.false.
6388 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
6389 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
6390 ixpmin1=max(1,ixp1-3)
6391 ixpmax1=min(ixp1+3,numxi1)
6392 ixpmin2=max(1,ixp2-3)
6393 ixpmax2=min(ixp2+3,numxi2)
6394 do ixp1=ixpmin1,ixpmax1
6395 do ixp2=ixpmin2,ixpmax2
6396 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6397 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6398 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6399 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6400 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
6401 wlb(ixp1,ixp2,1)=wlb(ixp1,ixp2,1)+tbsubc*factor
6402 wlb(ixp1,ixp2,2)=wlb(ixp1,ixp2,2)+pbsubc*factor
6412 deallocate(nh_dummy)
6419 double precision,
intent(in) :: Rl
6420 double precision,
intent(inout) :: A,B,C,D
6422 double precision :: sinO,cosO,sinO2,cosO2,tmp
6427 coso=abs(dsqrt(coso2))
6428 tmp=log((1.d0+sino)/coso)
6430 b=-(1.d0-3.d0*sino2-(coso2/sino)*(1.d0+3.d0*sino2)*tmp)/8.d0
6431 c=4.d0/3.d0-coso-coso*coso2/3.d0
6432 d=(5.d0+sino2-(coso2/sino)*(5.d0-sino2)*tmp)/8.d0
6438 double precision,
intent(in) :: Rl,Rin,Ne,A,B,C,D
6439 double precision,
intent(inout) :: fluxTB,fluxPB
6441 double precision :: const,u,Bt,Br,PB,TB,sinchi2
6444 const=1.24878d-25/(1.d0-u/3.d0)
6446 bt=const*(c+u*(d-c))
6447 pb=const*sinchi2*((a+u*(b-a)))
6456 double precision,
intent(in) :: x_sph(1:3)
6457 double precision,
intent(inout) :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
6459 unitv_r(1)=dsin(x_sph(2))*dcos(x_sph(3))
6460 unitv_r(2)=dsin(x_sph(2))*dsin(x_sph(3))
6461 unitv_r(3)=dcos(x_sph(2))
6462 unitv_theta(1)=dcos(x_sph(2))*dcos(x_sph(3))
6463 unitv_theta(2)=dcos(x_sph(2))*dsin(x_sph(3))
6464 unitv_theta(3)=-dsin(x_sph(2))
6465 unitv_phi(1)=-dsin(x_sph(3))
6466 unitv_phi(2)=dcos(x_sph(3))
6471 subroutine output_data(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,datatype)
6475 integer,
intent(in) :: qunit,nXO1,nXO2,nWO
6476 double precision,
intent(in) :: dxO1(nxO1),dxO2(nxO2)
6477 double precision,
intent(in) :: xO1(nXO1),xO2(nxO2)
6478 double precision,
intent(inout) :: wO(nXO1,nXO2,nWO)
6479 character(20),
intent(in) :: datatype
6481 integer :: nPiece,nP1,nP2,nC1,nC2,nWC
6482 integer :: piece_nmax1,piece_nmax2,ix1,ix2,j,ipc,ixc1,ixc2
6483 double precision :: uniform_tol
6484 double precision,
allocatable :: xC(:,:,:,:),wC(:,:,:,:),dxC(:,:,:,:)
6491 if (abs(wo(ix1,ix2,j))<smalldouble) wo(ix1,ix2,j)=zero
6498 if (datatype==
'image_euv' .or. datatype==
'image_sxr')
then
6500 piece_nmax1=block_nx2
6501 piece_nmax2=block_nx3
6503 piece_nmax1=block_nx3
6504 piece_nmax2=block_nx1
6506 piece_nmax1=block_nx1
6507 piece_nmax2=block_nx2
6509 else if (datatype==
'spectrum_euv')
then
6512 piece_nmax2=block_nx1
6514 piece_nmax2=block_nx2
6516 piece_nmax2=block_nx3
6523 loopn1:
do j=piece_nmax1,1,-1
6524 if(mod(nxo1,j)==0)
then
6529 loopn2:
do j=piece_nmax2,1,-1
6530 if(mod(nxo2,j)==0)
then
6543 case(
'EIvtuCCmpi',
'ESvtuCCmpi',
'SIvtuCCmpi',
'WIvtuCCmpi')
6545 allocate(xc(npiece,nc1,nc2,2))
6546 allocate(dxc(npiece,nc1,nc2,2))
6547 allocate(wc(npiece,nc1,nc2,nwo))
6551 ix1=mod(ipc-1,np1)*nc1+ixc1
6552 ix2=floor(1.0*(ipc-1)/np1)*nc2+ixc2
6553 xc(ipc,ixc1,ixc2,1)=xo1(ix1)
6554 xc(ipc,ixc1,ixc2,2)=xo2(ix2)
6555 dxc(ipc,ixc1,ixc2,1)=dxo1(ix1)
6556 dxc(ipc,ixc1,ixc2,2)=dxo2(ix2)
6558 wc(ipc,ixc1,ixc2,j)=wo(ix1,ix2,j)
6565 deallocate(xc,dxc,wc)
6566 case(
'EIvtiCCmpi',
'ESvtiCCmpi',
'SIvtiCCmpi',
'WIvtiCCmpi')
6568 (maxval(abs(dxo1(:)-dxo1(1)))>uniform_tol*max(one,abs(dxo1(1))) .or. &
6569 maxval(abs(dxo2(:)-dxo2(1)))>uniform_tol*max(one,abs(dxo2(1)))))
then
6570 call mpistop(
"vti needs uniform dat-resolution image grids")
6572 call write_image_vticc(qunit,xo1,xo2,dxo1,dxo2,wo,nxo1,nxo2,nwo,nc1,nc2)
6575 call mpistop(
"Error in synthesize emission: Unknown convert_type")
6581 subroutine write_image_vticc(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,nC1,nC2)
6585 integer,
intent(in) :: qunit,nXO1,nXO2,nWO,nC1,nC2
6586 double precision,
intent(in) :: xO1(nXO1),xO2(nxO2)
6587 double precision,
intent(in) :: dxO1(nxO1),dxO2(nxO2)
6588 double precision,
intent(in) :: wO(nXO1,nXO2,nWO)
6590 double precision :: origin(1:3), spacing(1:3)
6591 integer :: wholeExtent(1:6)
6593 integer :: ixC1,ixC2
6597 character (70) :: subname,wname,vname,nameL,nameS
6598 character (len=std_len) :: filename
6599 logical :: sph_datres_no_doppler
6602 origin(1)=xo1(1)-0.5d0*dxo1(1)
6603 origin(2)=xo2(1)-0.5d0*dxo2(1)
6614 inquire(qunit,opened=fileopen)
6615 if(.not.fileopen)
then
6620 write(filename,
'(a,i4.4,a)') trim(
filename_euv),filenr,
".vti"
6622 write(filename,
'(a,i4.4,a)') trim(
filename_sxr),filenr,
".vti"
6628 open(qunit,file=filename,status=
'unknown',form=
'formatted')
6632 write(qunit,
'(a)')
'<?xml version="1.0"?>'
6633 write(qunit,
'(a)',advance=
'no')
'<VTKFile type="ImageData"'
6634 write(qunit,
'(a)')
' version="0.1" byte_order="LittleEndian">'
6635 write(qunit,
'(a,3(1pe14.6),a,6(i10),a,3(1pe14.6),a)')
' <ImageData Origin="',&
6636 origin,
'" WholeExtent="',wholeextent,
'" Spacing="',spacing,
'">'
6638 write(qunit,
'(a)')
'<FieldData>'
6639 write(qunit,
'(2a)')
'<DataArray type="Float32" Name="TIME" ',&
6640 'NumberOfTuples="1" format="ascii">'
6642 write(qunit,
'(a)')
'</DataArray>'
6643 write(qunit,
'(a)')
'</FieldData>'
6645 write(qunit,
'(a,6(i10),a)')
'<Piece Extent="',wholeextent,
'">'
6646 write(qunit,
'(a)')
'<CellData>'
6657 if (trim(
emission_model)==
'pseudo_current' .and. iw==1) vname=
'pseudo_current'
6658 if (trim(
emission_model)==
'radio_ff' .and. iw==1) vname=
'radio_brightness_temperature'
6660 if (iw==2 .and.
dat_resolution .and. (.not. sph_datres_no_doppler) .and. &
6666 ((
dat_resolution .and. ((sph_datres_no_doppler .and. iw==2) .or. &
6667 ((.not. sph_datres_no_doppler) .and. iw==3))) .or. &
6681 vname=
'absorption_fraction'
6692 if (iw==1)
write(vname,
'(a)')
'B'
6693 if (iw==2)
write(vname,
'(a)')
'pB'
6701 write(qunit,
'(a,a,a)')&
6702 '<DataArray type="Float64" Name="',trim(vname),
'" format="ascii">'
6703 write(qunit,
'(200(1pe14.6))') ((wo(ixc1,ixc2,iw),ixc1=1,nxo1),ixc2=1,nxo2)
6704 write(qunit,
'(a)')
'</DataArray>'
6706 write(qunit,
'(a)')
'</CellData>'
6707 write(qunit,
'(a)')
'</Piece>'
6709 write(qunit,
'(a)')
'</ImageData>'
6710 write(qunit,
'(a)')
'</VTKFile>'
6720 integer,
intent(in) :: qunit
6721 integer,
intent(in) :: nPiece,nC1,nC2,nWC
6722 double precision,
intent(in) :: xC(nPiece,nC1,nC2,2),dxC(nPiece,nc1,nc2,2)
6723 double precision,
intent(in) :: wC(nPiece,nC1,nC2,nWC)
6724 character(20),
intent(in) :: datatype
6727 double precision :: xP(nPiece,nC1+1,nC2+1,2)
6730 character (70) :: subname,wname,vname,nameL,nameS
6731 character (len=std_len) :: filename
6732 integer :: ixC1,ixC2,ixP,ix1,ix2,j
6733 integer :: nc,np,icel,VTK_type
6734 logical :: sph_datres_no_doppler
6745 if (ix1<np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1,1,1)-0.5d0*dxc(ixp,ix1,1,1)
6746 if (ix1==np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1-1,1,1)+0.5d0*dxc(ixp,ix1-1,1,1)
6747 if (ix2<np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2,2)-0.5d0*dxc(ixp,1,ix2,2)
6748 if (ix2==np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2-1,2)+0.5d0*dxc(ixp,1,ix2-1,2)
6753 inquire(qunit,opened=fileopen)
6754 if(.not.fileopen)
then
6758 if (datatype==
'image_euv')
then
6759 write(filename,
'(a,i4.4,a)') trim(
filename_euv),filenr,
".vtu"
6760 else if (datatype==
'image_sxr')
then
6761 write(filename,
'(a,i4.4,a)') trim(
filename_sxr),filenr,
".vtu"
6762 else if (datatype==
'image_whitelight')
then
6764 else if (datatype==
'spectrum_euv')
then
6767 open(qunit,file=filename,status=
'unknown',form=
'formatted')
6770 write(qunit,
'(a)')
'<?xml version="1.0"?>'
6771 write(qunit,
'(a)',advance=
'no')
'<VTKFile type="UnstructuredGrid"'
6772 write(qunit,
'(a)')
' version="0.1" byte_order="LittleEndian">'
6773 write(qunit,
'(a)')
'<UnstructuredGrid>'
6774 write(qunit,
'(a)')
'<FieldData>'
6775 write(qunit,
'(2a)')
'<DataArray type="Float32" Name="TIME" ',&
6776 'NumberOfTuples="1" format="ascii">'
6778 write(qunit,
'(a)')
'</DataArray>'
6779 write(qunit,
'(a)')
'</FieldData>'
6781 write(qunit,
'(a,i7,a,i7,a)') &
6782 '<Piece NumberOfPoints="',np,
'" NumberOfCells="',nc,
'">'
6783 write(qunit,
'(a)')
'<CellData>'
6785 if (datatype==
'image_euv')
then
6794 if (trim(
emission_model)==
'pseudo_current') vname=
'pseudo_current'
6795 if (trim(
emission_model)==
'radio_ff') vname=
'radio_brightness_temperature'
6798 if (j==2 .and.
dat_resolution .and. (.not. sph_datres_no_doppler) .and. &
6804 ((
dat_resolution .and. ((sph_datres_no_doppler .and. j==2) .or. &
6805 ((.not. sph_datres_no_doppler) .and. j==3))) .or. &
6819 vname=
'absorption_fraction'
6821 else if (datatype==
'image_sxr')
then
6829 else if (datatype==
'image_whitelight')
then
6830 write(vname,
'(a)')
'whitelight'
6831 else if (datatype==
'spectrum_euv')
then
6838 write(qunit,
'(a,a,a)')&
6839 '<DataArray type="Float64" Name="',trim(vname),
'" format="ascii">'
6840 write(qunit,
'(200(1pe14.6))') ((wc(ixp,ixc1,ixc2,j),ixc1=1,nc1),ixc2=1,nc2)
6841 write(qunit,
'(a)')
'</DataArray>'
6843 write(qunit,
'(a)')
'</CellData>'
6844 write(qunit,
'(a)')
'<Points>'
6845 write(qunit,
'(a)')
'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
6850 write(qunit,
'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
6852 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
6854 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
6858 write(qunit,
'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
6860 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
6862 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
6865 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
6869 write(qunit,
'(a)')
'</DataArray>'
6870 write(qunit,
'(a)')
'</Points>'
6872 write(qunit,
'(a)')
'<Cells>'
6873 write(qunit,
'(a)')
'<DataArray type="Int32" Name="connectivity" format="ascii">'
6876 write(qunit,
'(4(i7))') ix1-1+(ix2-1)*np1,ix1+(ix2-1)*np1,&
6877 ix1-1+ix2*np1,ix1+ix2*np1
6880 write(qunit,
'(a)')
'</DataArray>'
6882 write(qunit,
'(a)')
'<DataArray type="Int32" Name="offsets" format="ascii">'
6884 write(qunit,
'(i7)') icel*(2**2)
6886 write(qunit,
'(a)')
'</DataArray>'
6888 write(qunit,
'(a)')
'<DataArray type="Int32" Name="types" format="ascii">'
6892 write(qunit,
'(i2)') vtk_type
6894 write(qunit,
'(a)')
'</DataArray>'
6895 write(qunit,
'(a)')
'</Cells>'
6896 write(qunit,
'(a)')
'</Piece>'
6898 write(qunit,
'(a)')
'</UnstructuredGrid>'
6899 write(qunit,
'(a)')
'</VTKFile>'
6905 double precision,
intent(in) :: vec1(1:3),vec2(1:3)
6906 double precision,
intent(out) :: res
6908 res=vec1(1)*vec2(1)+vec1(2)*vec2(2)+vec1(3)*vec2(3)
6913 double precision,
intent(in) :: vec_in1(1:3),vec_in2(1:3)
6914 double precision,
intent(out) :: vec_out(1:3)
6916 vec_out(1)=vec_in1(2)*vec_in2(3)-vec_in1(3)*vec_in2(2)
6917 vec_out(2)=vec_in1(3)*vec_in2(1)-vec_in1(1)*vec_in2(3)
6918 vec_out(3)=vec_in1(1)*vec_in2(2)-vec_in1(2)*vec_in2(1)
6924 double precision :: LOS_psi
6925 double precision :: vec_car(1:3),vec_z(1:3),vec_temp1(1:3),vec_temp2(1:3)
6926 double precision :: vec_LOS_sph(1:3),vec_xI1_sph(1:3),vec_xI2_sph(1:3)
6930 vec_los(2)=dpi*los_theta/180.d0
6941 if (los_theta==zero)
then
6943 vec_temp1(2)=dpi/2.d0
6944 vec_temp1(3)=dpi*los_phi/180.d0
6958 los_psi=dpi*image_rotate/180.d0
6959 vec_xi1=vec_temp1*cos(los_psi)-vec_temp2*sin(los_psi)
6960 vec_xi2=vec_temp2*cos(los_psi)+vec_temp1*sin(los_psi)
6971 vec_los_sph(2:3)=vec_los_sph(2:3)*180.d0/dpi
6972 vec_xi1_sph(2:3)=vec_xi1_sph(2:3)*180.d0/dpi
6973 vec_xi2_sph(2:3)=vec_xi2_sph(2:3)*180.d0/dpi
6975 if (mype==0)
write(*,
'(a,f3.1,f6.1,f6.1,a)')
' ray direction (spherical): [',vec_los_sph(1),vec_los_sph(2),vec_los_sph(3),
']'
6976 if (mype==0)
write(*,
'(a,f3.1,f6.1,f6.1,a)')
' xI1 direction (spherical): [',vec_xi1_sph(1),vec_xi1_sph(2),vec_xi1_sph(3),
']'
6977 if (mype==0)
write(*,
'(a,f3.1,f6.1,f6.1,a)')
' xI2 direction (spherical): [',vec_xi2_sph(1),vec_xi2_sph(2),vec_xi2_sph(3),
']'
6983 double precision,
intent(in) :: vec_sph(1:3)
6984 double precision,
intent(inout) :: vec_car(1:3)
6986 vec_car(1)=vec_sph(1)*dsin(vec_sph(2))*dcos(vec_sph(3))
6987 vec_car(2)=vec_sph(1)*dsin(vec_sph(2))*dsin(vec_sph(3))
6988 vec_car(3)=vec_sph(1)*dcos(vec_sph(2))
6994 double precision,
intent(in) :: vec_car(1:3)
6995 double precision,
intent(inout) :: vec_sph(1:3)
6997 vec_sph(1)=dsqrt(vec_car(1)**2+vec_car(2)**2+vec_car(3)**2)
6998 vec_sph(2)=dacos(vec_car(3)/vec_sph(1))
6999 vec_sph(3)=atan2(vec_car(2),vec_car(1))
7005 double precision :: LOS_psi
7006 double precision :: vec_z(1:3),vec_temp1(1:3),vec_temp2(1:3)
7009 vec_los(1)=-cos(dpi*los_phi/180.d0)*sin(dpi*los_theta/180.d0)
7010 vec_los(2)=-sin(dpi*los_phi/180.d0)*sin(dpi*los_theta/180.d0)
7011 vec_los(3)=-cos(dpi*los_theta/180.d0)
7017 if (los_theta==zero)
then
7018 vec_xi1(1)=cos(dpi*los_phi/180.d0)
7019 vec_xi1(2)=sin(dpi*los_phi/180.d0)
7027 los_psi=dpi*image_rotate/180.d0
7028 vec_xi1=vec_temp1*cos(los_psi)-vec_temp2*sin(los_psi)
7029 vec_xi2=vec_temp2*cos(los_psi)+vec_temp1*sin(los_psi)
7043 double precision,
intent(in) :: x_3D_sph(1:3)
7044 double precision,
intent(inout) :: x_image(1:2)
7045 double precision :: res,res_origin
7046 double precision :: x_3D(1:3)
7057 double precision,
intent(in) :: x_3D(1:3)
7058 double precision,
intent(inout) :: x_image(1:2)
7059 double precision :: res,res_origin
7063 x_image(1)=res-res_origin
7066 x_image(2)=res-res_origin
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
universal constants as specified in cgs units
PI (partial-ionisation) ionisation-degree backend for the eos% family.
subroutine, public ionization_state_tp(t, p, rfactor, iz_h, iz_he)
Equation of state for AMRVAC, handled through a single eos_container object.
Module with geometry-related routines (e.g., divergence, curl)
integer, parameter spherical
integer, parameter cartesian
subroutine curlvector(qvec, ixil, ixol, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
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.
logical output_absorption_fraction
output absorption fraction for thick/thin EUV synthesis when available
double precision radio_beam_fwhm
Gaussian radio beam full width at half maximum in arcsec.
integer snapshotini
Resume from the snapshot with this index.
character(len=std_len) filename_euv
Base file name for synthetic EUV emission output.
logical instrument_postprocess
Post-process dat-resolution EUV images onto the instrument pixel grid.
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.
integer, parameter rpxmin
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 radio_frequency
Observing frequency for radio free-free synthesis in Hz.
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
double precision, dimension(:), allocatable, parameter d
logical slab
Cartesian geometry or not.
double precision radio_beam_pixel_size
Output pixel size for radio beam post-processing in arcsec; <=0 uses FWHM/3.
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
integer npe
The number of MPI tasks.
logical output_tau
output optical-depth map for synthetic emission when available
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.
integer radsyn_segment_batch_factor
Maximum ray segments per pixel batch, as a factor of radsyn_pixel_batch; <=0 uses memory budget....
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)
character(len=std_len) dat_resolution_mode
Data-resolution image spacing: nominal or minimum actual cell size.
character(len=std_len) radiation_transfer
Synthetic emission transfer mode: thin or thick.
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
integer, parameter rpxmax
integer radsyn_pixel_batch
Number of image pixels processed in one ray-segment MPI batch.
logical radsyn_verbose
Print synthetic-emission ray-tracing profiling counters.
logical big_image
big image
double precision instrument_resolution_factor
times for enhancing spatial resolution for EUV image/spectra
double precision radsyn_segment_memory_mb
Approximate per-rank temporary memory budget, in MiB, for automatic ray-segment batch sizing.
double precision spectrum_window_min
spectral window
integer refine_max_level
Maximal number of AMR levels.
character(len=std_len) ray_method
Synthetic emission ray traversal method.
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
character(len=std_len) emission_model
Synthetic emission physical model selector.
integer, dimension(:,:), allocatable node
integer radsyn_segment_comm_factor
Maximum ray segments per segmented MPI all-to-all round, as a factor of radsyn_pixel_batch.
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...
double precision, dimension(1:3) vec_los
subroutine get_goes_flux_grid(ixil, ixol, w, x, dv, xboxl, xbl, fl, eflux_grid)
subroutine get_minimum_datresol_spacing(dxi)
subroutine integrate_spectra_cartesian(igrid, wl, dwlg, xs, dxsg, spectra, numwl, numxs, fl)
subroutine sph_cart_to_coord(pos, sph)
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
subroutine normalize_euv_doppler(ni1, ni2, euv, dpl, unitv)
subroutine get_sph_intersection_image_bounds(ximin1, ximax1, ximin2, ximax2)
double precision, dimension(1:60) f_263
integer function sph_locate_index(value, faces, imin, imax)
subroutine get_euv_image(qunit, fl)
double precision, dimension(1:60) t_eis1
subroutine postprocess_radio_beam_image(nsrc1, nsrc2, xsrc1, xsrc2, dxsrc1, dxsrc2, bright, nout1, nout2, xout1, xout2, dxout1, dxout2, wout, numwout, tau, brightthin)
subroutine solve_euv_saha_charge_state(nh, te, rhe, ne_guess, x_hii, x_heii, x_heiii)
subroutine get_sxr(ixil, ixol, w, x, fl, flux, el, eu)
integer function radsyn_euv_num_outputs(has_doppler, has_thick)
subroutine collect_euv_cart_dda_segments(ixil, ixol, source, opacity, sourcev, pixel_id, ray_origin, xface1, xface2, xface3, t_enter, t_exit, segments, nseg, capacity)
subroutine get_unit_vector_spherical(x_sph, unitv_r, unitv_theta, unitv_phi)
subroutine collect_euv_sph_intersection_segments(ixil, ixol, source, opacity, pixel_id, ray_origin, ximg1, ximg2, rface, thetaface, phiface, rface2, theta_cos, phi_sin, phi_cos, segments, nseg, capacity)
double precision, dimension(1:101) f_131
recursive subroutine quicksort_segment_indices(segments, idx, ilo, ihi)
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
subroutine sph_block_pixel_range(rface, thetaface, phiface, ixol, nxi1, nxi2, xi1, xi2, dxi, ixpmin1, ixpmax1, ixpmin2, ixpmax2, has_pixels)
subroutine cart_dda_advance_axis(ray_origin_axis, ray_dir_axis, faces, imin, imax, idx, step, tmax, done)
subroutine get_pseudo_current(igrid, ixil, ixol, w, source)
logical function segment_is_valid(segments, is, nvars)
double precision function transfer_attenuation(tau)
double precision, dimension(1:101) t_aia
double precision function exp_clamped(argument)
subroutine write_image_vtucc(qunit, xc, wc, dxc, npiece, nc1, nc2, nwc, datatype)
logical function sph_segment_visible(pos, ximg1, ximg2)
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)
subroutine cart_dda_block_pixel_range(box_min, box_max, nxif1, nxif2, xif1, xif2, ixpmin1, ixpmax1, ixpmin2, ixpmax2, has_pixels)
subroutine sph_add_t_fixed(tvals, nt, t)
double precision, dimension(1:60) t_eis2
subroutine sph_add_theta_intersections(ray_origin, ray_dir, thetaface, tvals, nt, capacity)
subroutine integrate_euv_cart_dda_thick_datresol(nxif1, nxif2, xif1, xif2, fl, euv, dpl, tau, euvthin)
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)
subroutine acc_euv_cart_dda(ixil, ixol, source, sourcev, ray_origin, xface1, xface2, xface3, t_enter, t_exit, euvp, dplp)
subroutine sph_add_phi_intersection(ray_origin, ray_dir, phiface, tvals, nt, capacity)
double precision, dimension(1:3) vec_xi1
subroutine append_cart_dda_segment(segments, nseg, capacity, pixel_id, tseg, jds, kds, jvds)
subroutine radsyn_get_segment_batch_limits(pixel_batch_target, segment_batch_target, segment_comm_target)
subroutine get_spectrum(qunit, datatype, fl)
subroutine ray_box_intersection_cart(ray_origin, ray_dir, box_min, box_max, hit, t_enter, t_exit)
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)
logical function radsyn_euv_has_doppler_output()
integer function sph_locate_index_desc(value, faces, imin, imax)
subroutine init_vectors_spherical()
subroutine insertion_sort_segment_indices(segments, idx, ilo, ihi)
subroutine get_nominal_datresol_spacing(dxi)
subroutine write_image_vticc(qunit, xo1, xo2, dxo1, dxo2, wo, nxo1, nxo2, nwo, nc1, nc2)
subroutine get_sxr_image(qunit, fl)
subroutine check_synthetic_emission_options(datatype)
subroutine fill_euv_absorption_fraction(ni1, ni2, euv, euvthin, smallflux, absorption, cap_to_one)
subroutine sph_sort_unique_t(tvals, nt)
subroutine init_vectors_cartesian()
double precision, dimension(1:60) f_192
subroutine sort_segment_indices_near_to_far(segments, idx, nidx)
integer function cart_dda_locate_index(pos, faces, imin, imax)
subroutine pack_euv_image_outputs(ni1, ni2, euv, wi, smallflux, has_doppler, has_thick, dpl, tau, euvthin, cap_absorption)
integer function segment_pixel_owner(pixel_id)
subroutine get_whitelight_thomson(rl, rin, ne, a, b, c, d, fluxtb, fluxpb)
subroutine get_image_datresol(qunit, datatype, fl)
subroutine collect_euv_sph_dda_interval(ixil, ixol, source, opacity, pixel_id, ray_origin, ximg1, ximg2, rface2, theta_cos, phiface, phi_sin, phi_cos, t_enter, t_exit, segments, nseg, capacity, ok)
subroutine get_goes_sxr_flux(xboxl, fl, eflux)
subroutine get_native_datresol_spacing(dxi)
double precision function interpolate_response_value(temperature, t_table, f_table, n_table, log_temperature, log_response)
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)
character(len=std_len) ray_method_active
subroutine integrate_sxr_datresol(igrid, nxif1, nxif2, xif1, xif2, dxif1, dxif2, fl, sxr)
subroutine apply_temperature_response(ixil, ixol, te, flux, t_table, f_table, n_table, log_temperature, log_response)
double precision, dimension(1:101) f_335
subroutine sph_try_theta_exit_candidate(t, theta_face_cos, ray_origin, tnow, texit, epsray, tnext, found)
subroutine get_radio_ff_source_opacity(ixil, ixol, w, x, fl, source, kappa)
subroutine sph_next_cell_exit(ray_origin, rface2, theta_cos, phiface, phi_sin, phi_cos, ixol, ix1, ix2, ix3, tnow, texit, epsray, tnext, found)
subroutine integrate_euv_thick_datresol(nxif1, nxif2, fl, euv, dpl, tau, euvthin)
subroutine output_data(qunit, xo1, xo2, dxo1, dxo2, wo, nxo1, nxo2, nwo, datatype)
subroutine integrate_euv_sph_intersection_thick(numxi1, numxi2, xi1, xi2, dxi, fl, euv, tau, euvthin)
subroutine sph_try_phi_exit_candidate(t, phi_face_sin, phi_face_cos, ray_origin, tnow, texit, epsray, tnext, found)
double precision, dimension(1:41) t_iris
subroutine sph_locate_cell_fast(pos, rface2, theta_cos, phiface, ixol, ix1, ix2, ix3, inside)
subroutine get_euv_spectrum(qunit, fl)
double precision, dimension(1:101) f_211
subroutine build_sph_intersection_faces(ixil, ixol, x, dx, rface, thetaface, phiface)
subroutine collect_euv_sph_dda_segments(ixil, ixol, source, opacity, pixel_id, ray_origin, ximg1, ximg2, rface, thetaface, phiface, rface2, theta_cos, phi_sin, phi_cos, segments, nseg, capacity, fallback)
subroutine cart_dda_init_axis(ray_origin_axis, ray_dir_axis, faces, imin, imax, idx, step, tmax)
subroutine get_euv_hhe_opacity(wl, ixil, ixol, w, x, fl, kappa)
subroutine integrate_euv_sph_intersection_thin(numxi1, numxi2, xi1, xi2, dxi, fl, em)
subroutine integrate_euv_cart_dda_datresol(nxif1, nxif2, xif1, xif2, fl, euv, dpl)
subroutine get_euv_saha_fractions(te, ne, x_hii, x_heii, x_heiii)
subroutine get_whitelight_image(qunit, fl)
subroutine get_euv(wl, ixil, ixol, w, x, fl, flux)
subroutine sph_add_sphere_intersections(ray_origin, ray_dir, rface, tvals, nt, capacity)
subroutine build_cart_dda_faces(ixil, ixol, x, dx, xface1, xface2, xface3)
subroutine integrate_transfer_step_first_order(emissivity, opacity, path_length, intensity, tau)
double precision function pow10_clamped(exponent)
subroutine acc_euv_sph_intersection(ixil, ixol, source, ray_origin, ximg1, ximg2, rface, thetaface, phiface, euvp)
subroutine sph_locate_cell(pos, rface, thetaface, phiface, ixol, ix1, ix2, ix3, inside)
subroutine sph_try_exit_candidate(t, tnow, texit, epsray, tnext, found)
subroutine spherical_to_cartesian(vec_sph, vec_car)
subroutine sph_add_t(tvals, nt, capacity, t)
subroutine postprocess_euv_instrument_image(nsrc1, nsrc2, xsrc1, xsrc2, dxsrc1, dxsrc2, euv, dpl, nout1, nout2, xout1, xout2, dxout1, dxout2, wout, numwout, tau, euvthin)