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(out):: val1(ixI^S), val2(ixI^S)
392 procedure(
get_subr1),
pointer,
nopass :: get_rho => null()
393 procedure(
get_subr1),
pointer,
nopass :: get_pthermal => null()
394 procedure(
get_subr1),
pointer,
nopass :: get_var_rfactor => null()
405 double precision,
allocatable :: source(:^d&)
406 double precision,
allocatable :: opacity(:^d&)
407 double precision,
allocatable :: sourcev(:^d&)
408 double precision,
allocatable :: xface1(:),xface2(:),xface3(:)
409 double precision,
allocatable :: rface(:),thetaface(:),phiface(:)
410 double precision,
allocatable :: rface2(:),theta_cos(:),phi_sin(:),phi_cos(:)
411 double precision :: box_min(1:3)=0.d0
412 double precision :: box_max(1:3)=0.d0
417 logical :: has_pixels=.false.
429 character(len=*),
intent(in) :: datatype
432 call mpistop(
"bad radiation_transfer")
437 call mpistop(
"dat_resolution_mode must be nominal or minimum")
453 case(
'cart',
'cart_dda')
455 case(
'spherical',
'sph_intersection')
457 case(
'sph_dda',
'spherical_dda')
467 call mpistop(
"bad emission_model")
471 call mpistop(
"tau and absorption-fraction output need thick transfer")
475 call mpistop(
"radsyn_pixel_batch must be positive")
478 call mpistop(
"radsyn_segment_batch_factor must be non-negative")
481 call mpistop(
"radsyn_segment_memory_mb must be positive")
484 call mpistop(
"radsyn_segment_comm_factor must be positive")
489 call mpistop(
"instrument_postprocess currently needs dat-resolution EUV images")
492 call mpistop(
"instrument_postprocess is not yet supported for spherical rays")
495 call mpistop(
"instrument_postprocess currently supports only EUV AIA or radio_ff images")
498 call mpistop(
"radio_ff instrument_postprocess needs radio_beam_fwhm > 0 arcsec")
506 if (datatype /=
'image_euv' .and. datatype /=
'spectrum_euv')
then
507 call mpistop(
"emission_model=euv_aia is only valid for EUV synthesis")
510 if (datatype /=
'image_whitelight')
then
511 call mpistop(
"emission_model=white_light is only valid for white-light synthesis")
514 if (datatype /=
'image_euv')
then
515 call mpistop(
"emission_model=radio_ff currently reuses EUV-image convert types")
518 call mpistop(
"emission_model=radio_ff needs radio_frequency > 0")
520 case(
'pseudo_current')
521 if (datatype /=
'image_euv')
then
522 call mpistop(
"emission_model=pseudo_current is only valid for EUV-image convert types")
525 call mpistop(
"emission_model=pseudo_current currently supports only thin transfer")
530 if (datatype /=
'image_euv' .or. .not.
slab)
then
531 call mpistop(
"ray_method=cart needs Cartesian EUV slab images")
536 call mpistop(
"ray_method=spherical currently needs 3D spherical grids")
539 call mpistop(
"ray_method=spherical currently needs 3D spherical grids")
544 call mpistop(
"bad ray_method=spherical mode")
547 call mpistop(
"ray_method=spherical currently supports only EUV AIA emission")
549 if (xprobmin2<=1.
d-10 .or. xprobmax2>=dpi-1.
d-10)
then
550 call mpistop(
"ray_method=spherical does not support polar-axis crossing domains")
552 if (xprobmax3<=xprobmin3 .or. xprobmax3-xprobmin3>=2.d0*dpi-1.
d-10)
then
553 call mpistop(
"ray_method=spherical does not support phi-wrapping domains")
558 if (datatype /=
'image_euv')
then
559 call mpistop(
"thick transfer is only defined for EUV images")
564 if (.not.
slab)
call mpistop(
"cartesian thick EUV currently needs slab output")
566 call mpistop(
"thick EUV currently needs Cartesian dat_resolution output")
572 call mpistop(
"thick EUV currently needs x/y/z-aligned LOS")
579 double precision,
intent(in) :: emissivity,opacity,path_length
580 double precision,
intent(inout) :: intensity,tau
582 double precision :: dtau
584 if (path_length<=zero)
return
586 dtau=max(zero,opacity)*path_length
592 trim(emission_model)/=
'radio_ff' .and. &
597 logical,
intent(in) :: has_doppler,has_thick
600 if (has_doppler) num_outputs=num_outputs+1
601 if (has_thick .and. output_tau) num_outputs=num_outputs+1
602 if (has_thick .and. output_absorption_fraction) num_outputs=num_outputs+1
606 integer,
intent(in) :: nI1,nI2
607 double precision,
intent(in) :: EUV(nI1,nI2),unitv
608 double precision,
intent(inout) :: Dpl(nI1,nI2)
614 if (euv(ix1,ix2)/=zero)
then
615 dpl(ix1,ix2)=(dpl(ix1,ix2)/euv(ix1,ix2))*unitv
619 if (abs(dpl(ix1,ix2))<smalldouble) dpl(ix1,ix2)=zero
625 integer,
intent(in) :: nI1,nI2
626 double precision,
intent(in) :: EUV(nI1,nI2),EUVthin(nI1,nI2),smallflux
627 double precision,
intent(out) :: Absorption(nI1,nI2)
628 logical,
intent(in),
optional :: cap_to_one
631 logical :: cap_absorption
634 cap_absorption=.false.
635 if (
present(cap_to_one)) cap_absorption=cap_to_one
638 if (euvthin(ix1,ix2)>smallflux)
then
639 absorption(ix1,ix2)=max(zero,(euvthin(ix1,ix2)-euv(ix1,ix2))/euvthin(ix1,ix2))
640 if (cap_absorption) absorption(ix1,ix2)=min(one,absorption(ix1,ix2))
646 subroutine pack_euv_image_outputs(nI1,nI2,EUV,wI,smallflux,has_doppler,has_thick,Dpl,Tau,EUVthin,&
648 integer,
intent(in) :: nI1,nI2
649 double precision,
intent(in) :: EUV(nI1,nI2),smallflux
650 double precision,
intent(inout) :: wI(:,:,:)
651 logical,
intent(in) :: has_doppler,has_thick
652 double precision,
intent(in),
optional :: Dpl(nI1,nI2),Tau(nI1,nI2),EUVthin(nI1,nI2)
653 logical,
intent(in),
optional :: cap_absorption
656 double precision,
allocatable :: Absorption(:,:)
661 if (has_doppler)
then
662 if (.not.
present(dpl))
call mpistop(
"Doppler output requested without Doppler image")
666 if (has_thick .and. output_tau)
then
667 if (.not.
present(tau))
call mpistop(
"tau output requested without tau image")
671 if (has_thick .and. output_absorption_fraction)
then
672 if (.not.
present(euvthin))
call mpistop(
"absorption output requested without thin image")
673 allocate(absorption(ni1,ni2))
676 wi(:,:,iw)=absorption(:,:)
677 deallocate(absorption)
682 integer,
intent(out) :: pixel_batch_target,segment_batch_target,segment_comm_target
684 pixel_batch_target=max(1,radsyn_pixel_batch)
685 if (radsyn_segment_batch_factor>0)
then
686 segment_batch_target=max(128,radsyn_segment_batch_factor*pixel_batch_target)
688 segment_batch_target=max(128,int(min(dble(huge(segment_batch_target)),&
689 max(128.d0,radsyn_segment_memory_mb*1048576.d0/256.d0))))
691 segment_comm_target=max(128,radsyn_segment_comm_factor*pixel_batch_target)
695 double precision,
intent(in) :: tau
705 double precision,
intent(in) :: argument
707 if (argument<-700.d0)
then
709 else if (argument>700.d0)
then
717 double precision,
intent(in) :: exponent
719 if (exponent>300.d0)
then
721 else if (exponent<-300.d0)
then
729 double precision,
intent(in) :: temperature
730 integer,
intent(in) :: n_table
731 double precision,
intent(in) :: t_table(n_table),f_table(n_table)
732 logical,
intent(in) :: log_temperature,log_response
734 integer :: ilo,ihi,imid
735 double precision :: temp_lookup,response_lookup,flo,fhi
738 if (temperature<=zero)
return
739 if (log_temperature)
then
740 temp_lookup=log10(temperature)
742 temp_lookup=temperature
744 if (temp_lookup<t_table(1) .or. temp_lookup>t_table(n_table))
return
745 if (temp_lookup==t_table(n_table))
then
746 if (log_response)
then
747 response_lookup=log10(max(f_table(n_table),1.d-99))
749 response_lookup=f_table(n_table)
756 if (temp_lookup>=t_table(imid))
then
762 if (log_response)
then
763 flo=log10(max(f_table(ilo),1.d-99))
764 fhi=log10(max(f_table(ilo+1),1.d-99))
769 response_lookup=flo*(temp_lookup-t_table(ilo+1))/(t_table(ilo)-t_table(ilo+1))+&
770 fhi*(temp_lookup-t_table(ilo))/(t_table(ilo+1)-t_table(ilo))
773 if (log_response)
then
782 integer,
intent(in) :: ixI^L, ixO^L, n_table
783 double precision,
intent(in) :: Te(ixI^S),t_table(n_table),f_table(n_table)
784 double precision,
intent(inout) :: flux(ixI^S)
785 logical,
intent(in) :: log_temperature,log_response
788 double precision :: GT
790 {
do ix^db=ixomin^db,ixomax^db\}
792 flux(ix^d)=flux(ix^d)*gt
793 if (flux(ix^d)<zero) flux(ix^d)=zero
803 double precision,
intent(in) :: Te,Ne
804 double precision,
intent(out) :: x_HII,x_HeII,x_HeIII
806 double precision :: Pe,log_H21,log_He21,log_He32,log_He321
807 double precision :: logScaleHe,w_H21,term0,term1,term2,denHe
808 double precision,
parameter :: Xe_H21=13.6d0
809 double precision,
parameter :: Xe_He21=24.587d0
810 double precision,
parameter :: Xe_He32=54.416d0
820 log_h21=2.5d0*log10(te)-5040.d0*xe_h21/te-log10(pe)-0.48d0
821 log_he21=log10(4.d0)+2.5d0*log10(te)-5040.d0*xe_he21/te-log10(pe)-0.48d0
822 log_he32=2.5d0*log10(te)-5040.d0*xe_he32/te-log10(pe)-0.48d0
825 x_hii=w_h21/(1.d0+w_h21)
829 log_he321=log_he21+log_he32
830 logscalehe=max(
zero,log_he21,log_he321)
834 denhe=term0+term1+term2
847 double precision,
intent(in) :: nH,Te,rHe,Ne_guess
848 double precision,
intent(out) :: x_HII,x_HeII,x_HeIII
850 integer,
parameter :: max_iter=32
852 double precision :: Ne,Ne_lo,Ne_hi,Ne_new,residual,derivative
853 double precision :: e_He,de_HII_dNe,de_He_dNe
858 if (nh<=zero .or. te<=zero)
return
860 ne_lo=max(1.d-30*nh,1.d-100)
861 ne_hi=(1.d0+2.d0*rhe)*nh
862 ne=min(max(ne_guess,ne_lo),ne_hi)
866 e_he=x_heii+2.d0*x_heiii
867 residual=ne/nh-x_hii-rhe*e_he
868 if (abs(residual)<1.d-10)
exit
870 if (residual>zero)
then
878 de_hii_dne=-x_hii*(1.d0-x_hii)/ne
879 de_he_dne=(x_heii*(e_he-1.d0) &
880 +2.d0*x_heiii*(e_he-2.d0))/ne
881 derivative=1.d0/nh-de_hii_dne-rhe*de_he_dne
882 ne_new=ne-residual/derivative
883 if (.not.(ne_new>ne_lo .and. ne_new<ne_hi))
then
884 ne_new=0.5d0*(ne_lo+ne_hi)
898 integer,
intent(in) :: wl
899 integer,
intent(in) :: ixI^L, ixO^L
900 double precision,
intent(in) :: x(ixI^S,1:ndim)
901 double precision,
intent(in) :: w(ixI^S,1:nw)
903 double precision,
intent(out) :: kappa(ixI^S)
906 double precision :: pth(ixI^S),rho(ixI^S),Rfactor(ixI^S),Te(ixI^S)
907 double precision :: Ne(ixI^S),nH(ixI^S)
908 double precision :: wave_ratio,s_H1,s_He1,s_He2
909 double precision :: x_HII,x_HeII,x_HeIII,iz_H,iz_He,Rdummy
910 double precision :: N_H1,N_He1,N_He2
911 double precision,
parameter :: rHe_opacity=0.1d0
912 double precision,
parameter :: sigma_H1=5.16d-20, sigma_he1=9.25d-19, sigma_he2=7.17d-19
914 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
915 call fl%get_rho(w,x,ixi^l,ixo^l,rho)
916 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,rfactor)
918 {
do ix^db=ixomin^db,ixomax^db\}
919 if (rho(ix^d)>zero .and. rfactor(ix^d)>zero)
then
920 te(ix^d)=pth(ix^d)/(rho(ix^d)*rfactor(ix^d))*unit_temperature
927 call eos%get_ne_nH(ixi^l, ixo^l, w, ne, nh)
929 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
930 nh(ixo^s)=nh(ixo^s)*unit_numberdensity/1.d6
932 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
933 nh(ixo^s)=nh(ixo^s)*unit_numberdensity
936 wave_ratio=dble(wl)/171.d0
940 if (wl<=912) s_h1=wave_ratio**3*sigma_h1
941 if (wl<=504) s_he1=wave_ratio**2*sigma_he1
942 if (wl<=228) s_he2=wave_ratio**2.75d0*sigma_he2
945 {
do ix^db=ixomin^db,ixomax^db\}
946 if (te(ix^d)>zero .and. nh(ix^d)>zero)
then
947 select case (trim(eos%eos_type))
955 x_heii=iz_he*(1.d0-iz_he)
962 if (ne(ix^d)>zero)
then
964 x_hii,x_heii,x_heiii)
965 if (trim(eos%method)==
'analytic')
then
968 x_hii=min(one,max(zero,ne(ix^d)/nh(ix^d)))
970 x_hii=min(one,max(zero,ne(ix^d)/nh(ix^d) &
971 -eos%He_abundance*(x_heii+2.d0*x_heiii)))
975 rhe_opacity,ne(ix^d),x_hii,x_heii,x_heiii)
982 rhe_opacity,ne(ix^d),x_hii,x_heii,x_heiii)
985 n_h1=nh(ix^d)*(1.d0-x_hii)
986 n_he1=rhe_opacity*nh(ix^d)*(1.d0-x_heii-x_heiii)
987 n_he2=rhe_opacity*nh(ix^d)*x_heii
988 kappa(ix^d)=max(zero,n_h1*s_h1+n_he1*s_he1+n_he2*s_he2)
994 integer,
intent(in) :: igrid
995 integer,
intent(in) :: ixI^L, ixO^L
996 double precision,
intent(in) :: w(ixI^S,1:nw)
997 double precision,
intent(out) :: source(ixI^S)
999 integer :: ix^D,idir,idirmin,idirmin0
1000 double precision :: current(ixI^S,7-2*ndir:3)
1002 if (.not.
allocated(iw_mag))
then
1003 call mpistop(
"emission_model=pseudo_current needs magnetic-field variables")
1008 call curlvector(w(ixi^s,iw_mag(1:ndir)),ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
1010 current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+ps(igrid)%J0(ixo^s,idirmin0:3)
1014 {
do ix^db=ixomin^db,ixomax^db\}
1016 source(ix^d)=source(ix^d)+current(ix^d,idir)**2
1024 integer,
intent(in) :: ixI^L, ixO^L
1025 double precision,
intent(in) :: x(ixI^S,1:ndim)
1026 double precision,
intent(in) :: w(ixI^S,1:nw)
1028 double precision,
intent(out) :: source(ixI^S),kappa(ixI^S)
1031 double precision :: pth(ixI^S),Te(ixI^S),Ne(ixI^S)
1032 double precision :: nH_dummy(ixI^S),gff
1034 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
1035 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
1036 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
1037 te(ixo^s)=pth(ixo^s)/(ne(ixo^s)*te(ixo^s))*unit_temperature
1038 call eos%get_ne_nH(ixi^l,ixo^l,w,ne,nh_dummy)
1040 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
1042 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
1047 {
do ix^db=ixomin^db,ixomax^db\}
1048 if (te(ix^d)>zero .and. ne(ix^d)>zero)
then
1049 if (te(ix^d)<2.d5)
then
1050 gff=18.2d0+1.5d0*log(te(ix^d))-log(radio_frequency)
1052 gff=24.5d0+log(te(ix^d))-log(radio_frequency)
1055 kappa(ix^d)=9.78d-3*ne(ix^d)**2*gff/(radio_frequency**2*te(ix^d)**1.5d0)
1056 source(ix^d)=te(ix^d)*kappa(ix^d)
1061 subroutine get_line_info(wl,ion,mass,logTe,line_center,spatial_px,spectral_px,sigma_PSF,width_slit)
1073 integer,
intent(in) :: wl
1074 integer,
intent(out) :: mass
1075 character(len=30),
intent(out) :: ion
1076 double precision,
intent(out) :: logTe,line_center,spatial_px,spectral_px
1077 double precision,
intent(out) :: sigma_PSF,width_slit
1147 line_center=1354.1d0
1149 spectral_px=12.98
d-3
1156 line_center=262.976d0
1165 line_center=263.765d0
1174 line_center=192.028d0
1183 line_center=255.113d0
1189 call mpistop(
"No information about this line")
1203 integer,
intent(in) :: wl
1204 integer,
intent(in) :: ixI^L, ixO^L
1205 double precision,
intent(in) :: x(ixI^S,1:ndim)
1206 double precision,
intent(in) :: w(ixI^S,1:nw)
1208 double precision,
intent(out) :: flux(ixI^S)
1211 double precision :: pth(ixI^S),rho(ixI^S),Rfactor(ixI^S),Te(ixI^S)
1212 double precision :: Ne(ixI^S),nH(ixI^S)
1214 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
1215 call fl%get_rho(w,x,ixi^l,ixo^l,rho)
1216 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,rfactor)
1222 call eos%get_ne_nH(ixi^l, ixo^l, w, ne, nh)
1230 flux(ixo^s)=ne(ixo^s)*nh(ixo^s)
1258 call mpistop(
"Unknown wavelength")
1270 integer,
intent(in) :: ixI^L,ixO^L
1271 integer,
intent(in) :: El,Eu
1272 double precision,
intent(in) :: x(ixI^S,1:ndim)
1273 double precision,
intent(in) :: w(ixI^S,nw)
1275 double precision,
intent(out) :: flux(ixI^S)
1277 integer :: ix^D,ixO^D
1279 double precision :: I0,kb,keV,dE,Ei
1280 double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
1281 double precision :: Ne(ixI^S),gff(ixI^S),fi(ixI^S)
1282 double precision :: EM(ixI^S)
1288 nume=floor((eu-el)/de)
1289 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
1290 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
1291 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
1295 double precision :: nH_dummy(ixI^S)
1296 call eos%get_ne_nH(ixi^l, ixo^l, w, ne, nh_dummy)
1299 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
1300 em(ixo^s)=(ne(ixo^s))**2*1.d6
1302 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
1303 em(ixo^s)=(ne(ixo^s))**2
1305 kbt(ixo^s)=kb*te(ixo^s)/kev
1310 {
do ix^db=ixomin^db,ixomax^db\}
1311 if (kbt(ix^d)>0.01*ei)
then
1312 if(kbt(ix^d)<ei) gff(ix^d)=(kbt(ix^d)/ei)**0.4
1313 fi(ix^d)=(em(ix^d)*gff(ix^d))*
exp_clamped(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
1318 flux(ixo^s)=flux(ixo^s)+fi(ixo^s)*de
1320 flux(ixo^s)=flux(ixo^s)*i0
1327 double precision,
intent(in) :: xbox^L
1329 double precision,
intent(out) :: eflux
1331 double precision :: dxb^D,xb^L
1332 integer :: iigrid,igrid,j
1333 integer :: ixO^L,ixI^L,ix^D
1334 double precision :: eflux_grid,eflux_pe
1336 ^d&iximin^d=
ixglo^d;
1337 ^d&iximax^d=
ixghi^d;
1338 ^d&ixomin^d=ixmlo^d;
1339 ^d&ixomax^d=ixmhi^d;
1341 do iigrid=1,igridstail; igrid=igrids(iigrid);
1345 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)
1346 eflux_pe=eflux_pe+eflux_grid
1348 call mpi_allreduce(eflux_pe,eflux,1,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
1355 integer,
intent(in) :: ixI^L,ixO^L
1356 double precision,
intent(in) :: x(ixI^S,1:ndim),dV(ixI^S)
1357 double precision,
intent(in) :: w(ixI^S,nw)
1358 double precision,
intent(in) :: xbox^L,xb^L
1360 double precision,
intent(out) :: eflux_grid
1362 integer :: ix^D,ixO^D,ixb^L
1363 integer :: iE,numE,j,inbox
1364 double precision :: I0,kb,keV,dE,Ei,El,Eu,A_cgs
1365 double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
1366 double precision :: Ne(ixI^S),EM(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)
1395 double precision :: nH_dummy(ixI^S)
1396 call eos%get_ne_nH(ixi^l, ixb^l, w, ne, nh_dummy)
1399 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
1400 em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*(unit_length*1.d2)**3
1402 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
1403 em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*unit_length**3
1405 kbt(ixb^s)=kb*te(ixb^s)/kev
1410 {
do ix^db=ixbmin^db,ixbmax^db\}
1411 if (kbt(ix^d)>1.d-2*ei)
then
1412 if(kbt(ix^d)<ei)
then
1413 gff=(kbt(ix^d)/ei)**0.4
1417 fi=(em(ix^d)*gff)*
exp_clamped(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
1418 eflux_grid=eflux_grid+fi*de*ei
1422 eflux_grid=eflux_grid*kev*erg_si
1431 integer,
intent(in) :: qunit
1433 character(20) :: datatype
1436 character (30) :: ion
1437 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1438 double precision :: xslit,arcsec
1440 datatype=
'spectrum_euv'
1445 if (
mype==0) print *,
'###################################################'
1448 if (
mype==0) print *,
'Systhesizing EUV spectrum (observed by IRIS).'
1449 case (263,264,192,255)
1450 if (
mype==0) print *,
'Systhesizing EUV spectrum (observed by Hinode/EIS).'
1452 call mpistop(
'Wrong wavelength!')
1456 call mpistop(
'Wrong spectrum window!')
1459 if (
mype==0)
write(*,
'(a,f8.3,a)')
' Wavelength: ',linecent,
' Angstrom'
1460 if (
mype==0) print *,
'Unit of EUV flux: DN s^-1 pixel^-1'
1464 write(*,
'(a,f5.3,a,f5.1,a)')
' Supposed pixel: ',wlrsl,
' Angstrom x ',spacersl*725.0,
' km'
1465 print *,
'Unit of wavelength: Angstrom (0.1 nm) '
1467 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
1469 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
1471 write(*,
'(a,f8.1,a)')
' Supposed width of slit: ',wslit*725.0,
' km'
1476 print *,
'Unit of wavelength: Angstrom (0.1 nm) '
1478 write(*,
'(a,f5.3,a,f5.1,a)')
' Pixel: ',wlrsl,
' Angstrom x ',spacersl*725.0,
' km'
1479 print *,
'Unit of length: arcsec (~725 km)'
1480 write(*,
'(a,f8.1,a)')
' Location of slit: xI1 = ',
location_slit,
' arcsec'
1481 write(*,
'(a,f8.1,a)')
' Width of slit: ',wslit,
' arcsec'
1484 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
1486 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
1488 write(*,
'(a,f8.1,a)')
' Location of slit: xI1 = ',
location_slit,
' Unit_length'
1489 write(*,
'(a,f8.1,a)')
' Width of slit: ',wslit*725.d0,
' km'
1492 if (
mype==0) print *,
'Direction of the slit: parallel to xI2 vector'
1493 if (coordinate==cartesian .or. coordinate==spherical)
then
1496 call mpistop(
"EUV spectrum synthesis: support for sperical coordinates is to be added!")
1500 if (
mype==0) print *,
'###################################################'
1506 integer,
intent(in) :: qunit
1507 character(20),
intent(in) :: datatype
1510 integer :: numWL,numXS,iwL,ixS,numWI,numS
1511 double precision :: dwLg,xSmin,xSmax,wLmin,wLmax
1512 double precision,
allocatable :: wL(:),xS(:),dwL(:),dxS(:)
1513 double precision,
allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
1514 integer :: strtype,nstrb,nbb,nuni,nstr,bnx
1515 double precision :: qs,dxfirst,dxmid,lenstr
1517 integer :: iigrid,igrid,j,dir_loc
1518 double precision :: xbmin(1:ndim),xbmax(1:ndim)
1521 numwl=4*int((spectrum_window_max-spectrum_window_min)/(4.d0*dwlg))
1522 wlmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwlg*numwl/2
1523 wlmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwlg*numwl/2
1524 allocate(wl(numwl),dwl(numwl))
1527 wl(iwl)=wlmin+iwl*dwlg-half*dwlg
1530 select case(direction_slit)
1532 numxs=domain_nx1*2**(refine_max_level-1)
1537 strtype=stretch_type(1)
1538 nstrb=nstretchedblocks_baselevel(1)
1539 qs=qstretch_baselevel(1)
1540 if (mype==0) print *,
'Direction of the slit: x'
1542 numxs=domain_nx2*2**(refine_max_level-1)
1547 strtype=stretch_type(2)
1548 nstrb=nstretchedblocks_baselevel(2)
1549 qs=qstretch_baselevel(2)
1550 if (mype==0) print *,
'Direction of the slit: y'
1552 numxs=domain_nx3*2**(refine_max_level-1)
1557 strtype=stretch_type(3)
1558 nstrb=nstretchedblocks_baselevel(3)
1559 qs=qstretch_baselevel(3)
1560 if (mype==0) print *,
'Direction of the slit: z'
1562 call mpistop(
'Wrong direction_slit')
1565 allocate(xs(numxs),dxs(numxs),spectra(numwl,numxs),spectra_rc(numwl,numxs))
1567 allocate(wi(numwl,numxs,numwi))
1569 select case(strtype)
1571 dxs(:)=(xsmax-xsmin)/numxs
1573 xs(ixs)=xsmin+dxs(ixs)*(ixs-half)
1576 qs=qs**(one/2**(refine_max_level-1))
1577 dxfirst=(xsmax-xsmin)*(one-qs)/(one-qs**numxs)
1580 dxs(ixs)=dxfirst*qs**(ixs-1)
1581 xs(ixs)=dxs(1)/(one-qs)*(one-qs**(ixs-1))+half*dxs(ixs)
1587 lenstr=(xsmax-xsmin)/(2.d0+nuni*(one-qs)/(one-qs**nstr))
1588 dxfirst=(xsmax-xsmin)/(dble(nuni)+2.d0/(one-qs)*(one-qs**nstr))
1591 nstr=nstr*2**(refine_max_level-1)
1592 nuni=nuni*2**(refine_max_level-1)
1593 qs=qs**(one/2**(refine_max_level-1))
1594 dxfirst=lenstr*(one-qs)/(one-qs**nstr)
1595 dxmid=dxmid/2**(refine_max_level-1)
1597 if(nuni .gt. 0)
then
1598 do ixs=nstr+1,nstr+nuni
1600 xs(ixs)=lenstr+(dble(ixs)-0.5d0-nstr)*dxs(ixs)+xsmin
1605 dxs(ixs)=dxfirst*qs**(nstr-ixs)
1606 xs(ixs)=xsmin+lenstr-dxs(ixs)*half-dxfirst*(one-qs**(nstr-ixs))/(one-qs)
1609 do ixs=nstr+nuni+1,numxs
1610 dxs(ixs)=dxfirst*qs**(ixs-nstr-nuni-1)
1611 xs(ixs)=xsmax-lenstr+dxs(ixs)*half+dxfirst*(one-qs**(ixs-nstr-nuni-1))/(one-qs)
1614 call mpistop(
"unknown stretch type")
1617 if (los_phi==0 .and. los_theta==90 .and. direction_slit==2)
then
1620 else if (los_phi==0 .and. los_theta==90 .and. direction_slit==3)
then
1623 else if (los_phi==90 .and. los_theta==90 .and. direction_slit==1)
then
1626 else if (los_phi==90 .and. los_theta==90 .and. direction_slit==3)
then
1629 else if (los_theta==0 .and. direction_slit==1)
then
1632 else if (los_theta==0 .and. direction_slit==2)
then
1636 call mpistop(
'Wrong combination of LOS and slit direction!')
1639 if (dir_loc==1)
then
1640 if (location_slit>xprobmax1 .or. location_slit<xprobmin1)
then
1641 call mpistop(
'Wrong value for location_slit!')
1643 if(mype==0)
write(*,
'(a,f8.1,a)')
' Location of slit: x = ',location_slit,
' Unit_length'
1644 else if (dir_loc==2)
then
1645 if (location_slit>xprobmax2 .or. location_slit<xprobmin2)
then
1646 call mpistop(
'Wrong value for location_slit!')
1648 if(mype==0)
write(*,
'(a,f8.1,a)')
' Location of slit: y = ',location_slit,
' Unit_length'
1650 if (location_slit>xprobmax3 .or. location_slit<xprobmin3)
then
1651 call mpistop(
'Wrong value for location_slit!')
1653 if(mype==0)
write(*,
'(a,f8.1,a)')
' Location of slit: z = ',location_slit,
' Unit_length'
1658 do iigrid=1,igridstail; igrid=igrids(iigrid);
1659 ^d&xbmin(^d)=rnode(rpxmin^d_,igrid);
1660 ^d&xbmax(^d)=rnode(rpxmax^d_,igrid);
1661 if (location_slit>=xbmin(dir_loc) .and. location_slit<xbmax(dir_loc))
then
1667 call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
1668 mpi_sum,icomm,ierrmpi)
1671 if (spectra_rc(iwl,ixs)>smalldouble)
then
1672 wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
1679 call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
1681 deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
1688 integer,
intent(in) :: igrid,numWL,numXS,dir_loc
1690 double precision,
intent(in) :: wL(numWL),dwL(numWL)
1691 double precision,
intent(inout) :: spectra(numWL,numXS)
1693 integer :: direction_LOS
1694 integer :: ixO^L,ixI^L,ix^D,ixOnew
1695 double precision,
allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
1696 double precision :: wlc,wlwd
1699 double precision :: logTe,lineCent
1700 character (30) :: ion
1701 double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
1703 integer :: levelg,rft,ixSmin,ixSmax,iwL
1704 double precision :: flux_pix,dL
1706 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1708 if (los_phi==0 .and. los_theta==90)
then
1710 else if (los_phi==90 .and. los_theta==90)
then
1716 ^d&ixomin^d=ixmlo^d\
1717 ^d&ixomax^d=ixmhi^d\
1718 ^d&iximin^d=ixglo^d\
1719 ^d&iximax^d=ixghi^d\
1720 allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
1723 if (dir_loc==1)
then
1724 do ix1=ixomin1,ixomax1
1725 if (location_slit>=(ps(igrid)%x(ix^d,1)-
half*ps(igrid)%dx(ix^d,1)) .and. &
1726 location_slit<(ps(igrid)%x(ix^d,1)+
half*ps(igrid)%dx(ix^d,1)))
then
1732 else if (dir_loc==2)
then
1733 do ix2=ixomin2,ixomax2
1734 if (location_slit>=(ps(igrid)%x(ix^d,2)-
half*ps(igrid)%dx(ix^d,2)) .and. &
1735 location_slit<(ps(igrid)%x(ix^d,2)+
half*ps(igrid)%dx(ix^d,2)))
then
1742 do ix3=ixomin3,ixomax3
1743 if (location_slit>=(ps(igrid)%x(ix^d,3)-
half*ps(igrid)%dx(ix^d,3)) .and. &
1744 location_slit<(ps(igrid)%x(ix^d,3)+
half*ps(igrid)%dx(ix^d,3)))
then
1752 call get_euv(spectrum_wl,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
1753 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
1754 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
1755 v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
1756 call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
1757 call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
1758 te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
1761 levelg=ps(igrid)%level
1762 rft=2**(refine_max_level-levelg)
1764 {
do ix^d=ixomin^d,ixomax^d\}
1767 wlc=linecent*(1.d0+v(ix^d)*unit_velocity*1.d2/
const_c)
1769 wlc=linecent*(1.d0+v(ix^d)*unit_velocity/
const_c)
1771 wlwd=sqrt(
kb_cgs*te(ix^d)*unit_temperature/(mass*
mp_cgs))
1774 select case(direction_slit)
1776 ixsmin=(block_nx1*(node(pig1_,igrid)-1)+(ix1-ixomin1))*rft+1
1777 ixsmax=(block_nx1*(node(pig1_,igrid)-1)+(ix1-ixomin1+1))*rft
1779 ixsmin=(block_nx2*(node(pig2_,igrid)-1)+(ix2-ixomin2))*rft+1
1780 ixsmax=(block_nx2*(node(pig2_,igrid)-1)+(ix2-ixomin2+1))*rft
1782 ixsmin=(block_nx3*(node(pig3_,igrid)-1)+(ix3-ixomin3))*rft+1
1783 ixsmax=(block_nx3*(node(pig3_,igrid)-1)+(ix3-ixomin3+1))*rft
1786 select case(direction_los)
1788 dl=ps(igrid)%dx(ix^d,1)*unit_length
1790 dl=ps(igrid)%dx(ix^d,2)*unit_length
1792 dl=ps(igrid)%dx(ix^d,3)*unit_length
1794 if (si_unit) dl=dl*1.d2
1797 flux_pix=flux(ix^d)*wlrsl*dl*
exp_clamped(-(wl(iwl)-wlc)**2/(2*wlwd**2))/(sqrt(2*
dpi)*wlwd)
1799 flux_pix=flux_pix*wslit/spacersl
1800 spectra(iwl,ixsmin:ixsmax)=spectra(iwl,ixsmin:ixsmax)+flux_pix
1806 deallocate(flux,v,pth,te,rho)
1812 integer,
intent(in) :: qunit
1813 character(20),
intent(in) :: datatype
1816 integer :: numWL,numXS,iwL,ixS,numWI,ix^D
1817 double precision :: dwLg,dxSg,xSmin,xSmax,xScent,wLmin,wLmax
1818 double precision,
allocatable :: wL(:),xS(:),dwL(:),dxS(:)
1819 double precision,
allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
1820 double precision :: vec_cor(1:3),xI_cor(1:2)
1821 double precision :: res,r_loc,r_max
1824 character (30) :: ion
1825 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1826 double precision :: unitv,arcsec,RHESSI_rsl,pixel
1827 integer :: iigrid,igrid,i,j,numS
1828 double precision :: xLmin,xLmax,xslit
1830 if (coordinate==spherical)
then
1838 if (coordinate==spherical)
then
1839 xsmin=-abs(xprobmax1)
1840 xsmax=abs(xprobmax1)
1843 if (ix1==1) vec_cor(1)=xprobmin1
1844 if (ix1==2) vec_cor(1)=xprobmax1
1846 if (ix2==1) vec_cor(2)=xprobmin2
1847 if (ix2==2) vec_cor(2)=xprobmax2
1849 if (ix3==1) vec_cor(3)=xprobmin3
1850 if (ix3==2) vec_cor(3)=xprobmax3
1852 r_loc=(vec_cor(1)-x_origin(1))**2
1853 r_loc=r_loc+(vec_cor(2)-x_origin(2))**2
1854 r_loc=r_loc+(vec_cor(3)-x_origin(3))**2
1856 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1859 r_max=max(r_max,r_loc)
1863 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1867 xsmin=min(xsmin,xi_cor(2))
1868 xsmax=max(xsmax,xi_cor(2))
1879 xscent=(xsmin+xsmax)/2.d0
1883 arcsec=7.25d5/unit_length
1885 arcsec=7.25d7/unit_length
1887 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1888 dxsg=spacersl*arcsec
1889 numxs=ceiling((xsmax-xscent)/dxsg)
1890 xsmin=xscent-numxs*dxsg
1891 xsmax=xscent+numxs*dxsg
1894 numwl=2*int((spectrum_window_max-spectrum_window_min)/(2.d0*dwlg))
1895 wlmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwlg*numwl/2
1896 wlmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwlg*numwl/2
1897 allocate(wl(numwl),dwl(numwl),xs(numxs),dxs(numxs))
1899 allocate(wi(numwl,numxs,numwi),spectra(numwl,numxs),spectra_rc(numwl,numxs))
1901 wl(iwl)=wlmin+iwl*dwlg-half*dwlg
1905 xs(ixs)=xsmin+dxsg*(ixs-half)
1911 do iigrid=1,igridstail; igrid=igrids(iigrid);
1913 if (ix1==1) vec_cor(1)=rnode(rpxmin1_,igrid)
1914 if (ix1==2) vec_cor(1)=rnode(rpxmax1_,igrid)
1916 if (ix2==1) vec_cor(2)=rnode(rpxmin2_,igrid)
1917 if (ix2==2) vec_cor(2)=rnode(rpxmax2_,igrid)
1919 if (ix3==1) vec_cor(3)=rnode(rpxmin3_,igrid)
1920 if (ix3==2) vec_cor(3)=rnode(rpxmax3_,igrid)
1922 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1926 xlmin=min(xlmin,xi_cor(1))
1927 xlmax=max(xlmax,xi_cor(1))
1933 if (activate_unit_arcsec)
then
1934 xslit=location_slit*arcsec
1938 if (xslit>=xlmin-wslit*arcsec .and. xslit<=xlmax+wslit*arcsec)
then
1944 call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
1945 mpi_sum,icomm,ierrmpi)
1948 if (spectra_rc(iwl,ixs)>smalldouble)
then
1949 wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
1956 if (activate_unit_arcsec)
then
1961 call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
1963 deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
1969 integer,
intent(in) :: igrid,numWL,numXS
1970 double precision,
intent(in) :: wL(numWL),xS(numXS)
1971 double precision,
intent(in) :: dwLg,dxSg
1972 double precision,
intent(inout) :: spectra(numWL,numXS)
1975 integer :: ixO^L,ixI^L,ix^D,ixOnew,j
1976 double precision,
allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
1977 double precision :: wlc,wlwd,res,dst_slit,xslit,arcsec
1978 double precision :: vloc(1:3),xloc(1:3),dxloc(1:3),xIloc(1:2),dxIloc(1:2)
1979 integer :: nSubC^D,iSubC^D,iwL,ixS,ixSmin,ixSmax,iwLmin,iwLmax,nwL
1980 double precision :: slit_width,dxSubC^D,xerf^L,fluxSubC
1981 double precision :: xSubC(1:3),xCent(1:2)
1984 double precision :: logTe,lineCent
1985 character (30) :: ion
1986 double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
1987 double precision :: sigma_wl,sigma_xs,factor
1990 arcsec=7.25d5/unit_length
1992 arcsec=7.25d7/unit_length
1994 if (activate_unit_arcsec)
then
1995 xslit=location_slit*arcsec
2000 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
2002 ^d&ixomin^d=ixmlo^d\
2003 ^d&ixomax^d=ixmhi^d\
2004 ^d&iximin^d=ixglo^d\
2005 ^d&iximax^d=ixghi^d\
2006 allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
2008 call get_euv(spectrum_wl,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
2009 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
2010 call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
2011 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
2012 call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
2013 te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
2014 {
do ix^d=ixomin^d,ixomax^d\}
2016 vloc(j)=ps(igrid)%w(ix^d,iw_mom(j))/rho(ix^d)
2024 slit_width=wslit*arcsec
2025 sigma_wl=sigma_psf*dwlg
2026 sigma_xs=sigma_psf*dxsg
2027 {
do ix^d=ixomin^d,ixomax^d\}
2028 if (flux(ix^d)>smalldouble)
then
2029 xloc(1:3)=ps(igrid)%x(ix^d,1:3)
2030 dxloc(1:3)=ps(igrid)%dx(ix^d,1:3)
2034 if (xiloc(1)>=xslit-half*(slit_width+dxiloc(1)) .and. &
2035 xiloc(1)<=xslit+half*(slit_width+dxiloc(1)))
then
2037 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi1(^d))/(slit_width/16.d0)));
2038 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi2(^d))/(dxsg/4.d0)));
2039 ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
2042 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length*1.d2/dxsg/dxsg
2043 wlc=linecent*(1.d0+v(ix^d)*unit_velocity*1.d2/const_c)
2045 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length/dxsg/dxsg
2046 wlc=linecent*(1.d0+v(ix^d)*unit_velocity/const_c)
2048 wlwd=sqrt(kb_cgs*te(ix^d)*unit_temperature/(mass*mp_cgs))
2049 wlwd=wlwd*linecent/const_c
2051 {
do isubc^d=1,nsubc^d\}
2052 ^d&xsubc(^d)=xloc(^d)-half*dxloc(^d)+(isubc^d-half)*dxsubc^d;
2054 dst_slit=abs(xcent(1)-xslit)
2055 if (dst_slit<=half*slit_width)
then
2056 ixs=floor((xcent(2)-(xs(1)-half*dxsg))/dxsg)+1
2058 ixsmax=min(ixs+3,numxs)
2059 iwl=floor((wlc-(wl(1)-half*dwlg))/dwlg)+1
2060 nwl=3*ceiling(wlwd/dwlg+1)
2061 iwlmin=max(1,iwl-nwl)
2062 iwlmax=min(iwl+nwl,numwl)
2064 do iwl=iwlmin,iwlmax
2065 do ixs=ixsmin,ixsmax
2066 xerfmin1=(wl(iwl)-half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
2067 xerfmax1=(wl(iwl)+half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
2068 xerfmin2=(xs(ixs)-half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
2069 xerfmax2=(xs(ixs)+half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
2070 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
2071 spectra(iwl,ixs)=spectra(iwl,ixs)+fluxsubc*factor
2081 deallocate(flux,v,pth,te)
2089 integer,
intent(in) :: qunit
2091 character(20) :: datatype
2094 character (30) :: ion
2095 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
2096 double precision :: t0,t1
2099 datatype=
'image_euv'
2104 print *,
'###################################################'
2105 print *,
'Systhesizing EUV image'
2106 write(*,
'(a,f8.3,a)')
' Wavelength: ',linecent,
' Angstrom'
2107 print *,
'Unit of EUV flux: DN s^-1 pixel^-1'
2112 call mpistop(
'EUV dat-resolution needs Cartesian or spherical native rays')
2114 print *,
'Data-resolution image requested; native output pixel sizes are reported below.'
2116 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2118 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2132 call mpistop(
'ERROR: Wrong LOS for synthesizing emission!')
2136 write(*,
'(a,f7.1,a,f7.1,a,f5.1,a,f5.1,a)')
' Pixel: ',spacersl*725.0,
' km x ',spacersl*725.0,
' km (', &
2137 spacersl,
' arcsec x ', spacersl,
' arcsec)'
2139 print *,
'Unit of length: arcsec (~725 km)'
2142 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2144 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2150 '] of the simulation box is located at [X=0,Y=0] of the image'
2152 else if (coordinate==spherical)
then
2153 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'
2156 call mpistop(
"EUV synthesis: this coordinate is not supported!")
2161 if (
mype==0) print *,
'time comsuming: ',t1-t0,
' s'
2162 if (
mype==0) print *,
'###################################################'
2169 integer,
intent(in) :: qunit
2171 character(20) :: datatype
2172 double precision :: RHESSI_rsl
2173 double precision :: t0,t1
2176 datatype=
'image_sxr'
2181 print *,
'###################################################'
2182 print *,
'Systhesizing SXR image (observed at 1 AU).'
2187 if (coordinate/=cartesian)
call mpistop(
'SXR synthesis: only cartesian is supported for .dat resolution!')
2189 print *,
'Unit of SXR flux: photons cm^-2 s^-1 pixel^-1'
2190 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, &
2191 ' Mm (', rhessi_rsl,
' arcsec x ', rhessi_rsl,
' arcsec)'
2193 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2195 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2205 call mpistop(
'ERROR: Wrong LOS for synthesizing emission!')
2209 print *,
'Unit of SXR flux: photons cm^-2 s^-1 pixel^-1'
2210 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, &
2211 ' Mm (', rhessi_rsl,
' arcsec x ', rhessi_rsl,
' arcsec)'
2213 print *,
'Unit of length: arcsec (~725 km)'
2216 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2218 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2222 if (coordinate==cartesian)
then
2224 '] of the simulation box is located at [X=0,Y=0] of the image'
2226 else if (coordinate==spherical)
then
2227 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'
2230 call mpistop(
"SXR synthesis: this coordinate is not supported!")
2235 if (
mype==0) print *,
'time comsuming:',t1-t0
2236 if (
mype==0) print *,
'###################################################'
2243 integer,
intent(in) :: qunit
2245 character(20) :: datatype
2246 double precision :: LASCO_rsl
2248 if (
mype==0) print *,
'###################################################'
2252 if (
mype==0) print *,
'Systhesizing white light image (observed by LASCO/C1).'
2255 if (
mype==0) print *,
'Systhesizing white light image (observed by LASCO/C2).'
2258 if (
mype==0) print *,
'Systhesizing white light image (observed by LASCO/C3).'
2260 call mpistop(
'Whitelight synthesis: instrument is not supported!')
2263 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 (', &
2264 lasco_rsl,
' arcsec x ', lasco_rsl,
' arcsec) '
2265 if (
mype==0) print *,
'Unit of white light flux: average Sun brightness'
2267 datatype=
'image_whitelight'
2272 print *,
'Unit of length: arcsec (~725 km)'
2275 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2277 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2282 if (coordinate==spherical)
then
2283 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'
2286 call mpistop(
"Whitelight synthesis: this coordinate is not supported!")
2289 if (
mype==0) print *,
'###################################################'
2294 EUV,Dpl,nOut1,nOut2,xOut1,xOut2,&
2295 dxOut1,dxOut2,wOut,numWOut,Tau,EUVthin)
2299 integer,
intent(in) :: nSrc1,nSrc2
2300 double precision,
intent(in) :: xSrc1(nSrc1),xSrc2(nSrc2)
2301 double precision,
intent(in) :: dxSrc1(nSrc1),dxSrc2(nSrc2)
2302 double precision,
intent(in) :: EUV(nSrc1,nSrc2),Dpl(nSrc1,nSrc2)
2303 integer,
intent(out) :: nOut1,nOut2,numWOut
2304 double precision,
allocatable,
intent(out) :: xOut1(:),xOut2(:),dxOut1(:),dxOut2(:)
2305 double precision,
allocatable,
intent(out) :: wOut(:,:,:)
2306 double precision,
intent(in),
optional :: Tau(nSrc1,nSrc2),EUVthin(nSrc1,nSrc2)
2308 integer :: mass,ixS1,ixS2,ixP1,ixP2,ixC1,ixC2,iw
2309 integer :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
2310 character(30) :: ion
2311 double precision :: logTe,lineCent,spaceRsl,wlRsl,sigma_PSF,wslit
2312 double precision :: arcsec,dxInst,xMin1,xMax1,xMin2,xMax2,xCent1,xCent2
2313 double precision :: sigma0,xerfmin1,xerfmax1,xerfmin2,xerfmax2
2314 double precision :: factor,weightSum,weightNorm,thinVal,tauVal
2315 double precision,
allocatable :: dplNum(:,:),thinOut(:,:),tauOut(:,:),tauWeight(:,:)
2323 dxinst=spacersl*arcsec
2324 if (dxinst<=
zero)
call mpistop(
"instrument_postprocess has non-positive pixel size")
2326 xmin1=minval(xsrc1-
half*dxsrc1)
2327 xmax1=maxval(xsrc1+
half*dxsrc1)
2328 xmin2=minval(xsrc2-
half*dxsrc2)
2329 xmax2=maxval(xsrc2+
half*dxsrc2)
2330 xcent1=
half*(xmin1+xmax1)
2331 xcent2=
half*(xmin2+xmax2)
2332 nout1=16*max(1,ceiling((xmax1-xmin1)/(16.d0*dxinst)))
2333 nout2=16*max(1,ceiling((xmax2-xmin2)/(16.d0*dxinst)))
2334 xmin1=xcent1-
half*dble(nout1)*dxinst
2335 xmin2=xcent2-
half*dble(nout2)*dxinst
2337 allocate(xout1(nout1),xout2(nout2),dxout1(nout1),dxout2(nout2))
2339 xout1(ixp1)=xmin1+dxinst*(dble(ixp1)-
half)
2343 xout2(ixp2)=xmin2+dxinst*(dble(ixp2)-
half)
2348 if (
present(tau) .and.
output_tau) numwout=numwout+1
2350 allocate(wout(nout1,nout2,numwout),dplnum(nout1,nout2))
2353 if (
present(euvthin))
then
2354 allocate(thinout(nout1,nout2))
2358 allocate(tauout(nout1,nout2),tauweight(nout1,nout2))
2363 sigma0=sigma_psf*dxinst
2368 if (
present(euvthin)) thinval=euvthin(ixs1,ixs2)
2369 if (
present(tau)) tauval=tau(ixs1,ixs2)
2373 ixc1=floor((xsrc1(ixs1)-(xout1(1)-
half*dxinst))/dxinst)+1
2374 ixc2=floor((xsrc2(ixs2)-(xout2(1)-
half*dxinst))/dxinst)+1
2375 ixpmin1=max(1,ixc1-3)
2376 ixpmax1=min(nout1,ixc1+3)
2377 ixpmin2=max(1,ixc2-3)
2378 ixpmax2=min(nout2,ixc2+3)
2381 do ixp1=ixpmin1,ixpmax1
2382 do ixp2=ixpmin2,ixpmax2
2383 xerfmin1=((xout1(ixp1)-
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2384 xerfmax1=((xout1(ixp1)+
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2385 xerfmin2=((xout2(ixp2)-
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2386 xerfmax2=((xout2(ixp2)+
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2387 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
2388 weightsum=weightsum+factor
2391 if (weightsum<=
zero) cycle
2393 do ixp1=ixpmin1,ixpmax1
2394 do ixp2=ixpmin2,ixpmax2
2395 xerfmin1=((xout1(ixp1)-
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2396 xerfmax1=((xout1(ixp1)+
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2397 xerfmin2=((xout2(ixp2)-
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2398 xerfmax2=((xout2(ixp2)+
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2399 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
2400 weightnorm=factor/weightsum
2401 wout(ixp1,ixp2,1)=wout(ixp1,ixp2,1)+euv(ixs1,ixs2)*weightnorm
2402 dplnum(ixp1,ixp2)=dplnum(ixp1,ixp2)+euv(ixs1,ixs2)*dpl(ixs1,ixs2)*weightnorm
2403 if (
present(euvthin)) thinout(ixp1,ixp2)=thinout(ixp1,ixp2)+thinval*weightnorm
2405 tauout(ixp1,ixp2)=tauout(ixp1,ixp2)+tauval*weightnorm
2406 tauweight(ixp1,ixp2)=tauweight(ixp1,ixp2)+weightnorm
2416 wout(ixp1,ixp2,2)=dplnum(ixp1,ixp2)/wout(ixp1,ixp2,1)
2418 wout(ixp1,ixp2,2)=
zero
2427 if (tauweight(ixp1,ixp2)>
zero)
then
2428 wout(ixp1,ixp2,iw)=tauout(ixp1,ixp2)/tauweight(ixp1,ixp2)
2430 wout(ixp1,ixp2,iw)=
zero
2440 wout(ixp1,ixp2,iw)=min(
one,max(
zero,(thinout(ixp1,ixp2)-wout(ixp1,ixp2,1))/thinout(ixp1,ixp2)))
2442 wout(ixp1,ixp2,iw)=
zero
2449 write(*,
'(a,2(i8,1x),a,2(i8,1x),a,1pe12.5)') &
2450 ' instrument_postprocess EUV grid src/out: ',nsrc1,nsrc2,
' -> ',nout1,nout2,
' dx=',dxinst
2454 if (
allocated(thinout))
deallocate(thinout)
2455 if (
allocated(tauout))
deallocate(tauout,tauweight)
2459 Bright,nOut1,nOut2,xOut1,xOut2,&
2460 dxOut1,dxOut2,wOut,numWOut,Tau,BrightThin)
2463 integer,
intent(in) :: nSrc1,nSrc2
2464 double precision,
intent(in) :: xSrc1(nSrc1),xSrc2(nSrc2)
2465 double precision,
intent(in) :: dxSrc1(nSrc1),dxSrc2(nSrc2)
2466 double precision,
intent(in) :: Bright(nSrc1,nSrc2)
2467 integer,
intent(out) :: nOut1,nOut2,numWOut
2468 double precision,
allocatable,
intent(out) :: xOut1(:),xOut2(:),dxOut1(:),dxOut2(:)
2469 double precision,
allocatable,
intent(out) :: wOut(:,:,:)
2470 double precision,
intent(in),
optional :: Tau(nSrc1,nSrc2),BrightThin(nSrc1,nSrc2)
2472 integer :: ixS1,ixS2,ixP1,ixP2,ixC1,ixC2,iw,nStencil
2473 integer :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
2474 double precision :: arcsec,beamPixel,beamSigma,xMin1,xMax1,xMin2,xMax2,xCent1,xCent2
2475 double precision :: distance1,distance2,weight,cellArea,thinVal,tauVal
2476 double precision,
allocatable :: norm(:,:),thinOut(:,:),tauOut(:,:),tauNorm(:,:)
2489 if (beamsigma<=zero .or. beampixel<=zero)
then
2490 call mpistop(
"radio beam postprocess has non-positive beam or pixel size")
2493 xmin1=minval(xsrc1-half*dxsrc1)
2494 xmax1=maxval(xsrc1+half*dxsrc1)
2495 xmin2=minval(xsrc2-half*dxsrc2)
2496 xmax2=maxval(xsrc2+half*dxsrc2)
2497 xcent1=half*(xmin1+xmax1)
2498 xcent2=half*(xmin2+xmax2)
2499 nout1=16*max(1,ceiling((xmax1-xmin1)/(16.d0*beampixel)))
2500 nout2=16*max(1,ceiling((xmax2-xmin2)/(16.d0*beampixel)))
2501 xmin1=xcent1-half*dble(nout1)*beampixel
2502 xmin2=xcent2-half*dble(nout2)*beampixel
2504 allocate(xout1(nout1),xout2(nout2),dxout1(nout1),dxout2(nout2))
2506 xout1(ixp1)=xmin1+beampixel*(dble(ixp1)-half)
2507 dxout1(ixp1)=beampixel
2510 xout2(ixp2)=xmin2+beampixel*(dble(ixp2)-half)
2511 dxout2(ixp2)=beampixel
2515 if (
present(tau) .and.
output_tau) numwout=numwout+1
2517 allocate(wout(nout1,nout2,numwout),norm(nout1,nout2))
2520 if (
present(brightthin))
then
2521 allocate(thinout(nout1,nout2))
2525 allocate(tauout(nout1,nout2),taunorm(nout1,nout2))
2530 nstencil=max(3,ceiling(4.d0*beamsigma/beampixel)+1)
2535 if (
present(brightthin)) thinval=brightthin(ixs1,ixs2)
2536 if (
present(tau)) tauval=tau(ixs1,ixs2)
2537 if (abs(bright(ixs1,ixs2))<=smalldouble .and. abs(thinval)<=smalldouble .and. &
2538 abs(tauval)<=smalldouble) cycle
2540 ixc1=floor((xsrc1(ixs1)-(xout1(1)-half*beampixel))/beampixel)+1
2541 ixc2=floor((xsrc2(ixs2)-(xout2(1)-half*beampixel))/beampixel)+1
2542 ixpmin1=max(1,ixc1-nstencil)
2543 ixpmax1=min(nout1,ixc1+nstencil)
2544 ixpmin2=max(1,ixc2-nstencil)
2545 ixpmax2=min(nout2,ixc2+nstencil)
2546 cellarea=max(smalldouble,dxsrc1(ixs1)*dxsrc2(ixs2))
2548 do ixp1=ixpmin1,ixpmax1
2549 distance1=xout1(ixp1)-xsrc1(ixs1)
2550 do ixp2=ixpmin2,ixpmax2
2551 distance2=xout2(ixp2)-xsrc2(ixs2)
2552 weight=
exp_clamped(-half*(distance1**2+distance2**2)/beamsigma**2)*cellarea
2553 wout(ixp1,ixp2,1)=wout(ixp1,ixp2,1)+bright(ixs1,ixs2)*weight
2554 norm(ixp1,ixp2)=norm(ixp1,ixp2)+weight
2555 if (
present(brightthin)) thinout(ixp1,ixp2)=thinout(ixp1,ixp2)+thinval*weight
2557 tauout(ixp1,ixp2)=tauout(ixp1,ixp2)+tauval*weight
2558 taunorm(ixp1,ixp2)=taunorm(ixp1,ixp2)+weight
2567 if (norm(ixp1,ixp2)>zero)
then
2568 wout(ixp1,ixp2,1)=wout(ixp1,ixp2,1)/norm(ixp1,ixp2)
2569 if (
present(brightthin)) thinout(ixp1,ixp2)=thinout(ixp1,ixp2)/norm(ixp1,ixp2)
2571 wout(ixp1,ixp2,1)=zero
2572 if (
present(brightthin)) thinout(ixp1,ixp2)=zero
2582 if (taunorm(ixp1,ixp2)>zero)
then
2583 wout(ixp1,ixp2,iw)=tauout(ixp1,ixp2)/taunorm(ixp1,ixp2)
2585 wout(ixp1,ixp2,iw)=zero
2594 if (thinout(ixp1,ixp2)>smalldouble)
then
2595 wout(ixp1,ixp2,iw)=min(one,max(zero,(thinout(ixp1,ixp2)-wout(ixp1,ixp2,1))/thinout(ixp1,ixp2)))
2597 wout(ixp1,ixp2,iw)=zero
2604 write(*,
'(a,2(i8,1x),a,2(i8,1x),a,2(1pe12.5,1x))') &
2605 ' radio_beam_postprocess grid src/out: ',nsrc1,nsrc2,
' -> ',nout1,nout2,&
2610 if (
allocated(thinout))
deallocate(thinout)
2611 if (
allocated(tauout))
deallocate(tauout,taunorm)
2620 integer,
intent(in) :: qunit
2621 character(20),
intent(in) :: datatype
2624 double precision :: dx^D
2625 integer :: numX^D,ix^D
2626 double precision,
allocatable :: EUV(:,:),EUVs(:,:),Dpl(:,:),Dpls(:,:)
2627 double precision,
allocatable :: EUVthin(:,:),Tau(:,:)
2628 double precision,
allocatable :: SXR(:,:),SXRs(:,:),wI(:,:,:)
2629 double precision,
allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:),dxIi
2630 integer :: numXI1,numXI2,numSI,numWI,iw
2631 double precision :: xI^L
2632 integer :: iigrid,igrid,i,j
2633 double precision,
allocatable :: xIF1(:),xIF2(:),dxIF1(:),dxIF2(:)
2634 double precision,
allocatable :: xIP1(:),xIP2(:),dxIP1(:),dxIP2(:),wIP(:,:,:)
2635 integer :: nXIF1,nXIF2
2636 integer :: nXIP1,nXIP2,numWIP
2637 double precision :: xIF^L
2638 double precision :: vec_cor(1:3),xI_cor(1:2),dxDDA,xIcent1,xIcent2
2640 double precision :: unitv,arcsec,RHESSI_rsl,length_to_km
2641 integer :: strtype^D,nstrb^D,nbb^D,nuni^D,nstr^D,bnx^D
2642 double precision :: qs^D,dxfirst^D,dxmid^D,lenstr^D
2643 logical :: has_doppler_output,has_thick_output
2656 xicent1=
half*(xifmin1+xifmax1)
2657 xicent2=
half*(xifmin2+xifmax2)
2658 nxif1=max(1,ceiling((xifmax1-xifmin1)/dxdda))
2659 nxif2=max(1,ceiling((xifmax2-xifmin2)/dxdda))
2660 xifmin1=xicent1-
half*dble(nxif1)*dxdda
2661 xifmax1=xicent1+
half*dble(nxif1)*dxdda
2662 xifmin2=xicent2-
half*dble(nxif2)*dxdda
2663 xifmax2=xicent2+
half*dble(nxif2)*dxdda
2674 if (
mype==0)
write(*,
'(a,a,a,1pe12.5,a,2(i8,1x))') &
2676 ' image-plane dx=',dxdda,
' n=',nxif1,nxif2
2679 if (ix1==1) vec_cor(1)=xprobmin1
2680 if (ix1==2) vec_cor(1)=xprobmax1
2682 if (ix2==1) vec_cor(2)=xprobmin2
2683 if (ix2==2) vec_cor(2)=xprobmax2
2685 if (ix3==1) vec_cor(3)=xprobmin3
2686 if (ix3==2) vec_cor(3)=xprobmax3
2688 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
2694 xifmin1=min(xifmin1,xi_cor(1))
2695 xifmax1=max(xifmax1,xi_cor(1))
2696 xifmin2=min(xifmin2,xi_cor(2))
2697 xifmax2=max(xifmax2,xi_cor(2))
2703 xicent1=
half*(xifmin1+xifmax1)
2704 xicent2=
half*(xifmin2+xifmax2)
2705 nxif1=max(1,ceiling((xifmax1-xifmin1)/dxdda))
2706 nxif2=max(1,ceiling((xifmax2-xifmin2)/dxdda))
2707 xifmin1=xicent1-
half*dble(nxif1)*dxdda
2708 xifmax1=xicent1+
half*dble(nxif1)*dxdda
2709 xifmin2=xicent2-
half*dble(nxif2)*dxdda
2710 xifmax2=xicent2+
half*dble(nxif2)*dxdda
2721 if (
mype==0)
write(*,
'(a,a,a,1pe12.5,a,2(i8,1x))') &
2723 ' image-plane dx=',dxdda,
' n=',nxif1,nxif2
2741 if (
mype==0)
write(*,
'(a)')
' LOS vector: [-1.00 0.00 0.00]'
2742 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [ 0.00 1.00 0.00]'
2743 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 0.00 1.00]'
2761 if (
mype==0)
write(*,
'(a)')
' LOS vector: [ 0.00 -1.00 0.00]'
2762 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [-1.00 0.00 0.00]'
2763 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 0.00 1.00]'
2781 if (
mype==0)
write(*,
'(a)')
' LOS vector: [ 0.00 0.00 -1.00]'
2782 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [ 1.00 0.00 0.00]'
2783 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 1.00 0.00]'
2785 allocate(xif1(nxif1),xif2(nxif2),dxif1(nxif1),dxif2(nxif2))
2788 select case(strtype1)
2790 dxif1(:)=(xifmax1-xifmin1)/nxif1
2792 xif1(ix1)=xifmin1+dxif1(ix1)*(ix1-
half)
2796 dxfirst1=(xifmax1-xifmin1)*(
one-qs1)/(
one-qs1**nxif1)
2799 dxif1(ix1)=dxfirst1*qs1**(ix1-1)
2800 xif1(ix1)=dxif1(1)/(
one-qs1)*(
one-qs1**(ix1-1))+
half*dxif1(ix1)
2805 nuni1=nbb1-nstrb1*bnx1
2806 lenstr1=(xifmax1-xifmin1)/(2.d0+nuni1*(
one-qs1)/(
one-qs1**nstr1))
2807 dxfirst1=(xifmax1-xifmin1)/(dble(nuni1)+2.d0/(
one-qs1)*(
one-qs1**nstr1))
2813 dxfirst1=lenstr1*(
one-qs1)/(
one-qs1**nstr1)
2816 if(nuni1 .gt. 0)
then
2817 do ix1=nstr1+1,nstr1+nuni1
2819 xif1(ix1)=lenstr1+(dble(ix1)-0.5d0-nstr1)*dxif1(ix1)+xifmin1
2824 dxif1(ix1)=dxfirst1*qs1**(nstr1-ix1)
2825 xif1(ix1)=xifmin1+lenstr1-dxif1(ix1)*
half-dxfirst1*(
one-qs1**(nstr1-ix1))/(
one-qs1)
2828 do ix1=nstr1+nuni1+1,nxif1
2829 dxif1(ix1)=dxfirst1*qs1**(ix1-nstr1-nuni1-1)
2830 xif1(ix1)=xifmax1-lenstr1+dxif1(ix1)*
half+dxfirst1*(
one-qs1**(ix1-nstr1-nuni1-1))/(
one-qs1)
2833 call mpistop(
"unknown stretch type")
2836 select case(strtype2)
2838 dxif2(:)=(xifmax2-xifmin2)/nxif2
2840 xif2(ix2)=xifmin2+dxif2(ix2)*(ix2-
half)
2844 dxfirst2=(xifmax2-xifmin2)*(
one-qs2)/(
one-qs2**nxif2)
2847 dxif2(ix2)=dxfirst2*qs2**(ix2-1)
2848 xif2(ix2)=dxif2(1)/(
one-qs2)*(
one-qs2**(ix2-1))+
half*dxif2(ix2)
2853 nuni2=nbb2-nstrb2*bnx2
2854 lenstr2=(xifmax2-xifmin2)/(2.d0+nuni2*(
one-qs2)/(
one-qs2**nstr2))
2855 dxfirst2=(xifmax2-xifmin2)/(dble(nuni2)+2.d0/(
one-qs2)*(
one-qs2**nstr2))
2861 dxfirst2=lenstr2*(
one-qs2)/(
one-qs2**nstr2)
2864 if(nuni2 .gt. 0)
then
2865 do ix2=nstr2+1,nstr2+nuni2
2867 xif2(ix2)=lenstr2+(dble(ix2)-0.5d0-nstr2)*dxif2(ix2)+xifmin2
2872 dxif2(ix2)=dxfirst2*qs2**(nstr2-ix2)
2873 xif2(ix2)=xifmin2+lenstr2-dxif2(ix2)*
half-dxfirst2*(
one-qs2**(nstr2-ix2))/(
one-qs2)
2876 do ix2=nstr2+nuni2+1,nxif2
2877 dxif2(ix2)=dxfirst2*qs2**(ix2-nstr2-nuni2-1)
2878 xif2(ix2)=xifmax2-lenstr2+dxif2(ix2)*
half+dxfirst2*(
one-qs2**(ix2-nstr2-nuni2-1))/(
one-qs2)
2881 call mpistop(
"unknown stretch type")
2884 if (
mype==0 .and. datatype==
'image_euv')
then
2892 write(*,
'(a,i8,a,i8)')
' Native data-resolution image grid: ',nxif1,
' x ',nxif2
2893 write(*,
'(a,f10.3,a,f10.3,a,f8.3,a,f8.3,a)') &
2894 ' Native xI1 pixel-size range: ',minval(dxif1)*length_to_km,
'--', &
2895 maxval(dxif1)*length_to_km,
' km (',minval(dxif1)/arcsec,
'--',maxval(dxif1)/arcsec,
' arcsec)'
2896 write(*,
'(a,f10.3,a,f10.3,a,f8.3,a,f8.3,a)') &
2897 ' Native xI2 pixel-size range: ',minval(dxif2)*length_to_km,
'--', &
2898 maxval(dxif2)*length_to_km,
' km (',minval(dxif2)/arcsec,
'--',maxval(dxif2)/arcsec,
' arcsec)'
2902 if (datatype==
'image_euv')
then
2911 allocate(wi(nxif1,nxif2,numwi))
2912 allocate(euv(nxif1,nxif2),dpl(nxif1,nxif2))
2914 allocate(euvthin(nxif1,nxif2),tau(nxif1,nxif2))
2923 if (has_doppler_output)
then
2928 allocate(euvs(nxif1,nxif2),dpls(nxif1,nxif2))
2937 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision, &
2942 do iigrid=1,igridstail; igrid=igrids(iigrid);
2946 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision, &
2948 call mpi_allreduce(dpls,dpl,numsi,mpi_double_precision, &
2951 if (has_doppler_output)
then
2955 deallocate(euvs,dpls)
2957 if (has_thick_output)
then
2958 if (has_doppler_output)
then
2960 has_thick_output,dpl=dpl,tau=tau,euvthin=euvthin)
2963 has_thick_output,tau=tau,euvthin=euvthin)
2965 else if (has_doppler_output)
then
2967 has_thick_output,dpl=dpl)
2977 euv,nxip1,nxip2,xip1,xip2,&
2978 dxip1,dxip2,wip,numwip,tau=tau,brightthin=euvthin)
2981 euv,nxip1,nxip2,xip1,xip2,&
2982 dxip1,dxip2,wip,numwip)
2986 euv,dpl,nxip1,nxip2,xip1,xip2,&
2987 dxip1,dxip2,wip,numwip,tau=tau,euvthin=euvthin)
2990 euv,dpl,nxip1,nxip2,xip1,xip2,&
2991 dxip1,dxip2,wip,numwip)
2993 call output_data(qunit,xip1,xip2,dxip1,dxip2,wip,nxip1,nxip2,numwip,datatype)
2994 deallocate(xip1,xip2,dxip1,dxip2,wip)
2996 call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
2999 deallocate(wi,euv,dpl,euvthin,tau)
3001 deallocate(wi,euv,dpl)
3006 if (datatype==
'image_sxr')
then
3014 allocate(wi(nxif1,nxif2,numwi))
3015 allocate(sxrs(nxif1,nxif2),sxr(nxif1,nxif2))
3018 do iigrid=1,igridstail; igrid=igrids(iigrid);
3022 call mpi_allreduce(sxrs,sxr,numsi,mpi_double_precision, &
3025 sxr=sxr*(rhessi_rsl*arcsec)**2
3033 call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
3034 deallocate(wi,sxr,sxrs)
3037 deallocate(xif1,xif2,dxif1,dxif2)
3044 integer,
intent(in) :: igrid,nXIF1,nXIF2
3045 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3046 double precision,
intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
3048 double precision,
intent(out) :: SXR(nXIF1,nXIF2)
3050 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
3051 double precision :: xb^L,xd^D
3052 double precision,
allocatable :: flux(:^D&),opacity(:^D&)
3053 double precision,
allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
3054 double precision,
allocatable :: SXRg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
3055 integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
3056 double precision :: SXRt,xc^L,xg^L,r2,area_1AU
3057 integer :: ixP^L,ixP^D
3058 integer :: direction_LOS
3068 ^d&ixomin^d=ixmlo^d\
3069 ^d&ixomax^d=ixmhi^d\
3070 ^d&iximin^d=
ixglo^d\
3071 ^d&iximax^d=
ixghi^d\
3075 allocate(flux(ixi^s))
3076 allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
3077 dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
3078 dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
3079 dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
3084 levelg=ps(igrid)%level
3088 select case(direction_los)
3099 allocate(sxrg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
3105 select case(direction_los)
3107 do ix2=ixomin2,ixomax2
3108 ixgmin1=(ix2-1)*rft+1
3110 do ix3=ixomin3,ixomax3
3111 ixgmin2=(ix3-1)*rft+1
3114 do ix1=ixomin1,ixomax1
3117 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
3121 do ix3=ixomin3,ixomax3
3122 ixgmin1=(ix3-1)*rft+1
3124 do ix1=ixomin1,ixomax1
3125 ixgmin2=(ix1-1)*rft+1
3128 do ix2=ixomin2,ixomax2
3131 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
3135 do ix1=ixomin1,ixomax1
3136 ixgmin1=(ix1-1)*rft+1
3138 do ix2=ixomin2,ixomax2
3139 ixgmin2=(ix2-1)*rft+1
3142 do ix3=ixomin3,ixomax3
3145 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
3155 select case(direction_los)
3157 ixgmin1=(ixomin2-1)*rft+1
3159 ixgmin2=(ixomin3-1)*rft+1
3162 ixgmin1=(ixomin3-1)*rft+1
3164 ixgmin2=(ixomin1-1)*rft+1
3167 ixgmin1=(ixomin1-1)*rft+1
3169 ixgmin2=(ixomin2-1)*rft+1
3173 select case(direction_los)
3175 ixpmin1=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3176 ixpmax1=
node(pig2_,igrid)*rft*block_nx2
3177 ixpmin2=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3178 ixpmax2=
node(pig3_,igrid)*rft*block_nx3
3180 ixpmin1=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3181 ixpmax1=
node(pig3_,igrid)*rft*block_nx3
3182 ixpmin2=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3183 ixpmax2=
node(pig1_,igrid)*rft*block_nx1
3185 ixpmin1=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3186 ixpmax1=
node(pig1_,igrid)*rft*block_nx1
3187 ixpmin2=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3188 ixpmax2=
node(pig2_,igrid)*rft*block_nx2
3190 xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
3191 xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
3192 dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
3193 dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
3194 sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
3195 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
3197 deallocate(flux,dxb1,dxb2,dxb3,sxrg,xg1,xg2,dxg1,dxg2)
3204 integer,
intent(in) :: igrid,nXIF1,nXIF2
3205 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3206 double precision,
intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
3208 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
3210 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
3211 double precision :: xb^L,xd^D
3212 double precision,
allocatable :: flux(:^D&),v(:^D&),rho(:^D&),opacity(:^D&)
3213 double precision,
allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
3214 double precision,
allocatable :: EUVg(:,:),Fvg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
3215 integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
3216 double precision :: EUVt,Fvt,xc^L,xg^L,r2
3217 integer :: ixP^L,ixP^D
3218 integer :: direction_LOS
3228 ^d&ixomin^d=ixmlo^d\
3229 ^d&ixomax^d=ixmhi^d\
3230 ^d&iximin^d=
ixglo^d\
3231 ^d&iximax^d=
ixghi^d\
3235 allocate(flux(ixi^s),v(ixi^s),rho(ixi^s),opacity(ixi^s))
3236 allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
3237 dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
3238 dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
3239 dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
3252 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
3253 v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
3258 levelg=ps(igrid)%level
3262 select case(direction_los)
3273 allocate(euvg(nxg1,nxg2),fvg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
3280 select case(direction_los)
3282 do ix2=ixomin2,ixomax2
3283 ixgmin1=(ix2-1)*rft+1
3285 do ix3=ixomin3,ixomax3
3286 ixgmin2=(ix3-1)*rft+1
3290 do ix1=ixomin1,ixomax1
3294 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
3295 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
3299 do ix3=ixomin3,ixomax3
3300 ixgmin1=(ix3-1)*rft+1
3302 do ix1=ixomin1,ixomax1
3303 ixgmin2=(ix1-1)*rft+1
3307 do ix2=ixomin2,ixomax2
3311 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
3312 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
3316 do ix1=ixomin1,ixomax1
3317 ixgmin1=(ix1-1)*rft+1
3319 do ix2=ixomin2,ixomax2
3320 ixgmin2=(ix2-1)*rft+1
3324 do ix3=ixomin3,ixomax3
3328 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
3329 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
3340 select case(direction_los)
3342 ixgmin1=(ixomin2-1)*rft+1
3344 ixgmin2=(ixomin3-1)*rft+1
3347 ixgmin1=(ixomin3-1)*rft+1
3349 ixgmin2=(ixomin1-1)*rft+1
3352 ixgmin1=(ixomin1-1)*rft+1
3354 ixgmin2=(ixomin2-1)*rft+1
3358 select case(direction_los)
3360 ixpmin1=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3361 ixpmax1=
node(pig2_,igrid)*rft*block_nx2
3362 ixpmin2=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3363 ixpmax2=
node(pig3_,igrid)*rft*block_nx3
3365 ixpmin1=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3366 ixpmax1=
node(pig3_,igrid)*rft*block_nx3
3367 ixpmin2=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3368 ixpmax2=
node(pig1_,igrid)*rft*block_nx1
3370 ixpmin1=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3371 ixpmax1=
node(pig1_,igrid)*rft*block_nx1
3372 ixpmin2=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3373 ixpmax2=
node(pig2_,igrid)*rft*block_nx2
3375 xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
3376 xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
3377 dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
3378 dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
3379 euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
3380 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
3381 dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
3382 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
3384 deallocate(flux,v,opacity,dxb1,dxb2,dxb3,euvg,fvg,xg1,xg2,dxg1,dxg2)
3393 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),box_min(1:3),box_max(1:3)
3394 logical,
intent(out) :: hit
3395 double precision,
intent(out) :: t_enter,t_exit
3398 double precision :: t1,t2,td
3404 if (abs(ray_dir(idir))<=smalldouble)
then
3405 if (ray_origin(idir)<box_min(idir) .or. ray_origin(idir)>box_max(idir))
then
3410 t1=(box_min(idir)-ray_origin(idir))/ray_dir(idir)
3411 t2=(box_max(idir)-ray_origin(idir))/ray_dir(idir)
3417 t_enter=max(t_enter,t1)
3418 t_exit=min(t_exit,t2)
3419 if (t_enter>=t_exit)
then
3428 integer,
intent(in) :: ixI^L, ixO^L
3429 double precision,
intent(in) :: x(ixI^S,1:ndim),dx(ixI^S,1:ndim)
3430 double precision,
allocatable,
intent(out) :: xface1(:),xface2(:),xface3(:)
3434 allocate(xface1(ixomin1:ixomax1+1),xface2(ixomin2:ixomax2+1),xface3(ixomin3:ixomax3+1))
3438 do ix1=ixomin1,ixomax1
3439 xface1(ix1)=x(ix^d,1)-half*dx(ix^d,1)
3441 xface1(ixomax1+1)=x(ixomax1,ixomin2,ixomin3,1)+half*dx(ixomax1,ixomin2,ixomin3,1)
3445 do ix2=ixomin2,ixomax2
3446 xface2(ix2)=x(ix^d,2)-half*dx(ix^d,2)
3448 xface2(ixomax2+1)=x(ixomin1,ixomax2,ixomin3,2)+half*dx(ixomin1,ixomax2,ixomin3,2)
3452 do ix3=ixomin3,ixomax3
3453 xface3(ix3)=x(ix^d,3)-half*dx(ix^d,3)
3455 xface3(ixomax3+1)=x(ixomin1,ixomin2,ixomax3,3)+half*dx(ixomin1,ixomin2,ixomax3,3)
3459 integer,
intent(in) :: imin,imax
3460 double precision,
intent(in) :: pos,faces(imin:imax+1)
3462 integer :: ilo,ihi,imid
3464 if (pos<=faces(imin))
then
3468 if (pos>=faces(imax+1))
then
3475 do while (ihi-ilo>1)
3477 if (pos>=faces(imid))
then
3483 idx=min(imax,max(imin,ilo))
3487 integer,
intent(in) :: imin,imax,idx
3488 double precision,
intent(in) :: ray_origin_axis,ray_dir_axis,faces(imin:imax+1)
3489 integer,
intent(out) :: step
3490 double precision,
intent(out) :: tMax
3492 if (ray_dir_axis>zero)
then
3494 tmax=(faces(idx+1)-ray_origin_axis)/ray_dir_axis
3495 else if (ray_dir_axis<zero)
then
3497 tmax=(faces(idx)-ray_origin_axis)/ray_dir_axis
3505 integer,
intent(in) :: imin,imax,step
3506 double precision,
intent(in) :: ray_origin_axis,ray_dir_axis,faces(imin:imax+1)
3507 integer,
intent(inout) :: idx
3508 double precision,
intent(inout) :: tMax
3509 logical,
intent(out) :: done
3513 if (idx<imin .or. idx>imax)
then
3518 tmax=(faces(idx+1)-ray_origin_axis)/ray_dir_axis
3519 else if (step<0)
then
3520 tmax=(faces(idx)-ray_origin_axis)/ray_dir_axis
3527 ray_origin,xface1,xface2,xface3,t_enter,t_exit,EUVp,Dplp)
3528 integer,
intent(in) :: ixI^L, ixO^L
3529 double precision,
intent(in) :: source(ixI^S),sourcev(ixI^S)
3530 double precision,
intent(in) :: ray_origin(1:3)
3531 double precision,
intent(in) :: xface1(ixOmin1:ixOmax1+1),xface2(ixOmin2:ixOmax2+1),&
3532 xface3(ixOmin3:ixOmax3+1)
3533 double precision,
intent(in) :: t_enter,t_exit
3534 double precision,
intent(inout) :: EUVp,Dplp
3536 integer :: ix^D,step(1:3)
3537 double precision :: pos(1:3),tMax(1:3),tNow,tNext,ds_cm,epsRay
3540 if (t_exit<=t_enter)
return
3541 epsray=max(1.d-12,1.d-10*abs(t_exit-t_enter))
3542 pos=ray_origin+(t_enter+epsray)*
vec_los
3553 tnext=min(t_exit,tmax(1),tmax(2),tmax(3))
3554 if (tnext>tnow)
then
3555 ds_cm=(tnext-tnow)*unit_length
3556 if (si_unit) ds_cm=ds_cm*1.d2
3557 euvp=euvp+source(ix^d)*ds_cm
3558 dplp=dplp+sourcev(ix^d)*ds_cm
3561 if (tnow>=t_exit-epsray)
exit
3563 if (tmax(1)<=tnow+epsray)
then
3567 if (tmax(2)<=tnow+epsray)
then
3571 if (tmax(3)<=tnow+epsray)
then
3579 double precision,
allocatable,
intent(inout) :: segments(:,:)
3580 integer,
intent(inout) :: nseg,capacity
3581 integer,
intent(in) :: pixel_id
3582 double precision,
intent(in) :: tseg,jds,kds,jvds
3584 double precision,
allocatable :: tmp(:,:)
3585 integer :: new_capacity
3587 if (capacity<=0)
then
3589 allocate(segments(5,capacity))
3590 else if (nseg>=capacity)
then
3591 new_capacity=2*capacity
3592 allocate(tmp(5,new_capacity))
3593 tmp(:,1:capacity)=segments(:,1:capacity)
3594 call move_alloc(tmp,segments)
3595 capacity=new_capacity
3599 segments(1,nseg)=dble(pixel_id)
3600 segments(2,nseg)=tseg
3601 segments(3,nseg)=jds
3602 segments(4,nseg)=kds
3603 segments(5,nseg)=jvds
3607 pixel_id,ray_origin,xface1,xface2,xface3,t_enter,t_exit,&
3608 segments,nseg,capacity)
3609 integer,
intent(in) :: ixI^L, ixO^L
3610 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S),sourcev(ixI^S)
3611 integer,
intent(in) :: pixel_id
3612 double precision,
intent(in) :: ray_origin(1:3)
3613 double precision,
intent(in) :: xface1(ixOmin1:ixOmax1+1),xface2(ixOmin2:ixOmax2+1),&
3614 xface3(ixOmin3:ixOmax3+1)
3615 double precision,
intent(in) :: t_enter,t_exit
3616 double precision,
allocatable,
intent(inout) :: segments(:,:)
3617 integer,
intent(inout) :: nseg,capacity
3619 integer :: ix^D,step(1:3)
3620 double precision :: pos(1:3),tMax(1:3),tNow,tNext,ds_cm,epsRay,tseg
3621 double precision :: jds,kds,jvds
3624 if (t_exit<=t_enter)
return
3625 epsray=max(1.d-12,1.d-10*abs(t_exit-t_enter))
3626 pos=ray_origin+(t_enter+epsray)*
vec_los
3637 tnext=min(t_exit,tmax(1),tmax(2),tmax(3))
3638 if (tnext>tnow)
then
3639 ds_cm=(tnext-tnow)*unit_length
3640 if (si_unit) ds_cm=ds_cm*1.d2
3641 jds=source(ix^d)*ds_cm
3642 kds=opacity(ix^d)*ds_cm
3643 jvds=sourcev(ix^d)*ds_cm
3644 if (jds/=zero .or. kds/=zero .or. jvds/=zero)
then
3645 tseg=half*(tnow+tnext)
3650 if (tnow>=t_exit-epsray)
exit
3652 if (tmax(1)<=tnow+epsray)
then
3656 if (tmax(2)<=tnow+epsray)
then
3660 if (tmax(3)<=tnow+epsray)
then
3668 double precision,
intent(in) :: segments(:,:)
3669 integer,
intent(inout) :: idx(:)
3670 integer,
intent(in) :: nidx
3681 double precision,
intent(in) :: segments(:,:)
3682 integer,
intent(inout) :: idx(:)
3683 integer,
intent(in) :: ilo,ihi
3690 do while (j>=ilo .and. segments(2,idx(j))>segments(2,key))
3699 double precision,
intent(in) :: segments(:,:)
3700 integer,
intent(inout) :: idx(:)
3701 integer,
intent(in) :: ilo,ihi
3704 double precision :: pivot
3706 if (ihi-ilo<=32)
then
3713 pivot=segments(2,idx((ilo+ihi)/2))
3715 do while (segments(2,idx(i))<pivot)
3718 do while (segments(2,idx(j))>pivot)
3736 integer,
intent(in) :: pixel_id
3738 owner=mod(pixel_id-1,npe)
3742 double precision,
intent(in) :: segments(:,:)
3743 integer,
intent(in) :: is,nvars
3749 if (segments(iv,is)/=segments(iv,is) .or. abs(segments(iv,is))>=1.d90)
then
3756 subroutine cart_dda_block_pixel_range(box_min,box_max,nXIF1,nXIF2,xIF1,xIF2,ixPmin1,ixPmax1,ixPmin2,ixPmax2,has_pixels)
3757 double precision,
intent(in) :: box_min(1:3),box_max(1:3)
3758 integer,
intent(in) :: nXIF1,nXIF2
3759 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3760 integer,
intent(out) :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
3761 logical,
intent(out) :: has_pixels
3764 double precision :: vec_cor(1:3),xI_cor(1:2)
3765 double precision :: xmin1,xmax1,xmin2,xmax2,dx1,dx2
3768 if (i1==1) vec_cor(1)=box_min(1)
3769 if (i1==2) vec_cor(1)=box_max(1)
3771 if (i2==1) vec_cor(2)=box_min(2)
3772 if (i2==2) vec_cor(2)=box_max(2)
3774 if (i3==1) vec_cor(3)=box_min(3)
3775 if (i3==2) vec_cor(3)=box_max(3)
3777 if (i1==1 .and. i2==1 .and. i3==1)
then
3783 xmin1=min(xmin1,xi_cor(1))
3784 xmax1=max(xmax1,xi_cor(1))
3785 xmin2=min(xmin2,xi_cor(2))
3786 xmax2=max(xmax2,xi_cor(2))
3793 dx1=abs(xif1(2)-xif1(1))
3795 dx1=max(one,abs(xmax1-xmin1))
3798 dx2=abs(xif2(2)-xif2(1))
3800 dx2=max(one,abs(xmax2-xmin2))
3803 ixpmin1=max(1,floor((xmin1-xif1(1))/dx1)+1-1)
3804 ixpmax1=min(nxif1,ceiling((xmax1-xif1(1))/dx1)+1+1)
3805 ixpmin2=max(1,floor((xmin2-xif2(1))/dx2)+1-1)
3806 ixpmax2=min(nxif2,ceiling((xmax2-xif2(1))/dx2)+1+1)
3807 has_pixels=ixpmin1<=ixpmax1 .and. ixpmin2<=ixpmax2
3813 integer,
intent(in) :: nXIF1,nXIF2
3814 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3816 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
3818 integer :: ixO^L,ixO^D,ixI^L,ix^D
3819 integer :: iigrid,igrid,ixP1,ixP2,numSI,ixPmin1,ixPmax1,ixPmin2,ixPmax2
3820 double precision :: box_min(1:3),box_max(1:3),ray_origin(1:3)
3821 double precision :: t_enter,t_exit,vlos
3822 double precision :: profile_local(2),profile_global(2)
3823 logical :: hit,has_pixels
3824 double precision,
allocatable :: source(:^D&),sourcev(:^D&),rho(:^D&),opacity(:^D&)
3825 double precision,
allocatable :: xface1(:),xface2(:),xface3(:)
3826 double precision,
allocatable :: EUVs(:,:),Dpls(:,:)
3828 allocate(euvs(nxif1,nxif2),dpls(nxif1,nxif2))
3833 do iigrid=1,igridstail; igrid=igrids(iigrid);
3834 ^d&ixomin^d=ixmlo^d\
3835 ^d&ixomax^d=ixmhi^d\
3836 ^d&iximin^d=
ixglo^d\
3837 ^d&iximax^d=
ixghi^d\
3839 box_min(1)=
rnode(rpxmin1_,igrid)
3840 box_min(2)=
rnode(rpxmin2_,igrid)
3841 box_min(3)=
rnode(rpxmin3_,igrid)
3842 box_max(1)=
rnode(rpxmax1_,igrid)
3843 box_max(2)=
rnode(rpxmax2_,igrid)
3844 box_max(3)=
rnode(rpxmax3_,igrid)
3847 ixpmin1,ixpmax1,ixpmin2,ixpmax2,has_pixels)
3848 if (.not. has_pixels)
then
3849 deallocate(xface1,xface2,xface3)
3853 allocate(source(ixi^s),sourcev(ixi^s),rho(ixi^s),opacity(ixi^s))
3863 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
3864 do ix1=ixomin1,ixomax1
3865 do ix2=ixomin2,ixomax2
3866 do ix3=ixomin3,ixomax3
3867 if (rho(ix^d)>smalldouble)
then
3868 vlos=(ps(igrid)%w(ix^d,iw_mom(1))*
vec_los(1)+&
3869 ps(igrid)%w(ix^d,iw_mom(2))*
vec_los(2)+&
3870 ps(igrid)%w(ix^d,iw_mom(3))*
vec_los(3))/rho(ix^d)
3871 sourcev(ix^d)=source(ix^d)*vlos
3877 deallocate(rho,opacity)
3879 do ixp1=ixpmin1,ixpmax1
3880 do ixp2=ixpmin2,ixpmax2
3882 profile_local(1)=profile_local(1)+one
3885 profile_local(2)=profile_local(2)+one
3887 ray_origin,xface1,xface2,xface3,t_enter,t_exit,euvs(ixp1,ixp2),dpls(ixp1,ixp2))
3892 deallocate(source,sourcev,xface1,xface2,xface3)
3896 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
3897 call mpi_allreduce(dpls,dpl,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
3898 call mpi_allreduce(profile_local,profile_global,2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
3900 write(*,
'(a,2(es12.5,1x))')
' cart_dda thin profile ray_tests ray_hits: ',profile_global
3902 deallocate(euvs,dpls)
3908 integer,
intent(in) :: nXIF1,nXIF2
3909 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3911 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
3912 double precision,
intent(out) :: Tau(nXIF1,nXIF2),EUVthin(nXIF1,nXIF2)
3914 integer,
parameter :: nSegVars=5
3915 integer :: ixO^L,ixO^D,ixI^L,ix^D
3916 integer :: iigrid,igrid,ixP1,ixP2,ipix,ipixStart,ipixEnd,nPixBatch,pixel_id
3917 integer :: nseg,capacity,totalCount,totalSeg,ipe,is,iseg,nidx,owner,isegDest,nsegBefore
3918 integer :: ixGlobal,iyGlobal,ixPmin1,ixPmax1,ixPmin2,ixPmax2,iFirst,iLast,iLocal
3919 integer :: nPixBatchTarget
3920 integer :: maxSegBatchTarget,maxSegCommTarget,maxNsegBatch,nPixTotal
3921 integer :: maxOwnerSegCount,maxOwnerSegCountLocal,segOffset,recvFill,totalRoundCount,totalRoundSeg
3922 integer,
allocatable :: sendCounts(:),recvCounts(:),sendDispls(:),recvDispls(:)
3923 integer,
allocatable :: roundSendCounts(:),roundRecvCounts(:)
3924 integer,
allocatable :: roundSendDispls(:),roundRecvDispls(:)
3925 integer,
allocatable :: ownerSegCounts(:),ownerOffsets(:),idx(:)
3926 integer,
allocatable :: bucketCounts(:),bucketOffsets(:),bucketFill(:)
3927 double precision :: ray_origin(1:3)
3928 double precision :: t_enter,t_exit,vlos,atten
3929 double precision :: profile_local(5),profile_global(5),profile_batch(5)
3930 logical :: hit,has_pixels,batchAccepted,batchReduced
3931 double precision,
allocatable :: rho(:^D&)
3932 double precision,
allocatable :: segments(:,:),segments_send(:,:),segments_recv(:,:)
3933 double precision,
allocatable :: segments_recv_round(:,:)
3934 double precision,
allocatable :: image_reduce(:,:)
3942 allocate(sendcounts(0:
npe-1),recvcounts(0:
npe-1),senddispls(0:
npe-1),recvdispls(0:
npe-1))
3943 allocate(roundsendcounts(0:
npe-1),roundrecvcounts(0:
npe-1))
3944 allocate(roundsenddispls(0:
npe-1),roundrecvdispls(0:
npe-1))
3945 allocate(ownersegcounts(0:
npe-1),owneroffsets(0:
npe-1))
3946 allocate(cache(igridstail))
3948 allocate(bucketcounts(npixbatchtarget),bucketoffsets(npixbatchtarget+1),&
3949 bucketfill(npixbatchtarget))
3951 do iigrid=1,igridstail; igrid=igrids(iigrid);
3952 ^d&ixomin^d=ixmlo^d\
3953 ^d&ixomax^d=ixmhi^d\
3954 ^d&iximin^d=
ixglo^d\
3955 ^d&iximax^d=
ixghi^d\
3957 cache(iigrid)%igrid=igrid
3958 allocate(cache(iigrid)%source(ixi^s),cache(iigrid)%opacity(ixi^s),&
3959 cache(iigrid)%sourcev(ixi^s),rho(ixi^s))
3960 cache(iigrid)%source=zero
3961 cache(iigrid)%opacity=zero
3962 cache(iigrid)%sourcev=zero
3965 cache(iigrid)%source,cache(iigrid)%opacity)
3967 call get_euv(
wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,cache(iigrid)%source)
3970 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
3971 do ix1=ixomin1,ixomax1
3972 do ix2=ixomin2,ixomax2
3973 do ix3=ixomin3,ixomax3
3974 if (rho(ix^d)>smalldouble)
then
3975 vlos=(ps(igrid)%w(ix^d,iw_mom(1))*
vec_los(1)+&
3976 ps(igrid)%w(ix^d,iw_mom(2))*
vec_los(2)+&
3977 ps(igrid)%w(ix^d,iw_mom(3))*
vec_los(3))/rho(ix^d)
3978 cache(iigrid)%sourcev(ix^d)=cache(iigrid)%source(ix^d)*vlos
3985 cache(iigrid)%box_min(1)=
rnode(rpxmin1_,igrid)
3986 cache(iigrid)%box_min(2)=
rnode(rpxmin2_,igrid)
3987 cache(iigrid)%box_min(3)=
rnode(rpxmin3_,igrid)
3988 cache(iigrid)%box_max(1)=
rnode(rpxmax1_,igrid)
3989 cache(iigrid)%box_max(2)=
rnode(rpxmax2_,igrid)
3990 cache(iigrid)%box_max(3)=
rnode(rpxmax3_,igrid)
3992 cache(iigrid)%xface1,cache(iigrid)%xface2,&
3993 cache(iigrid)%xface3)
3995 nxif1,nxif2,xif1,xif2,cache(iigrid)%ixPmin1,cache(iigrid)%ixPmax1,&
3996 cache(iigrid)%ixPmin2,cache(iigrid)%ixPmax2,cache(iigrid)%has_pixels)
3999 npixtotal=nxif1*nxif2
4001 do while (ipixstart<=npixtotal)
4002 ipixend=min(nxif1*nxif2,ipixstart+npixbatchtarget-1)
4003 npixbatch=ipixend-ipixstart+1
4004 batchaccepted=.false.
4005 batchreduced=.false.
4007 do while (.not. batchaccepted)
4012 do iigrid=1,igridstail; igrid=igrids(iigrid);
4013 ^d&ixomin^d=ixmlo^d\
4014 ^d&ixomax^d=ixmhi^d\
4015 ^d&iximin^d=
ixglo^d\
4016 ^d&iximax^d=
ixghi^d\
4018 ixpmin1=cache(iigrid)%ixPmin1
4019 ixpmax1=cache(iigrid)%ixPmax1
4020 ixpmin2=cache(iigrid)%ixPmin2
4021 ixpmax2=cache(iigrid)%ixPmax2
4022 has_pixels=cache(iigrid)%has_pixels
4023 if (.not. has_pixels) cycle
4025 do ixp2=ixpmin2,ixpmax2
4026 ifirst=max(ipixstart,(ixp2-1)*nxif1+ixpmin1)
4027 ilast=min(ipixend,(ixp2-1)*nxif1+ixpmax1)
4028 if (ifirst>ilast) cycle
4029 do ipix=ifirst,ilast
4030 ixp1=1+mod(ipix-1,nxif1)
4032 profile_batch(1)=profile_batch(1)+one
4034 cache(iigrid)%box_max,hit,t_enter,t_exit)
4036 profile_batch(2)=profile_batch(2)+one
4039 cache(iigrid)%opacity,cache(iigrid)%sourcev,&
4040 ipix,ray_origin,cache(iigrid)%xface1,&
4041 cache(iigrid)%xface2,cache(iigrid)%xface3,&
4043 segments,nseg,capacity)
4044 profile_batch(3)=profile_batch(3)+dble(nseg-nsegbefore)
4050 call mpi_allreduce(nseg,maxnsegbatch,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
4051 if (maxnsegbatch>maxsegbatchtarget .and. npixbatch>1)
then
4052 npixbatch=max(1,npixbatch/2)
4053 ipixend=ipixstart+npixbatch-1
4054 if (
allocated(segments))
deallocate(segments)
4057 batchaccepted=.true.
4061 profile_local=profile_local+profile_batch
4063 write(*,
'(a,3(i0,1x))')
' cart_dda thick adaptive batch: ',&
4064 ipixstart,ipixend,maxnsegbatch
4067 if (.not.
allocated(segments))
then
4069 allocate(segments(nsegvars,capacity))
4074 ownersegcounts(owner)=ownersegcounts(owner)+1
4076 sendcounts=nsegvars*ownersegcounts
4079 senddispls(ipe)=senddispls(ipe-1)+sendcounts(ipe-1)
4082 allocate(segments_send(nsegvars,max(1,nseg)))
4086 isegdest=senddispls(owner)/nsegvars+owneroffsets(owner)+1
4087 segments_send(:,isegdest)=segments(:,is)
4088 owneroffsets(owner)=owneroffsets(owner)+1
4091 call mpi_alltoall(sendcounts,1,mpi_integer,recvcounts,1,mpi_integer,
icomm,
ierrmpi)
4094 recvdispls(ipe)=recvdispls(ipe-1)+recvcounts(ipe-1)
4096 totalcount=sum(recvcounts)
4097 totalseg=totalcount/nsegvars
4098 profile_local(4)=profile_local(4)+dble(totalcount)
4099 allocate(segments_recv(nsegvars,max(1,totalseg)))
4102 maxownersegcountlocal=maxval(ownersegcounts)
4103 call mpi_allreduce(maxownersegcountlocal,maxownersegcount,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
4104 do segoffset=0,maxownersegcount-1,maxsegcommtarget
4106 roundsenddispls=senddispls
4108 if (ownersegcounts(ipe)>segoffset)
then
4109 roundsendcounts(ipe)=nsegvars*min(maxsegcommtarget,ownersegcounts(ipe)-segoffset)
4110 roundsenddispls(ipe)=senddispls(ipe)+nsegvars*segoffset
4114 call mpi_alltoall(roundsendcounts,1,mpi_integer,roundrecvcounts,1,mpi_integer,
icomm,
ierrmpi)
4115 roundrecvdispls(0)=0
4117 roundrecvdispls(ipe)=roundrecvdispls(ipe-1)+roundrecvcounts(ipe-1)
4119 totalroundcount=sum(roundrecvcounts)
4120 totalroundseg=totalroundcount/nsegvars
4121 allocate(segments_recv_round(nsegvars,max(1,totalroundseg)))
4123 call mpi_alltoallv(segments_send,roundsendcounts,roundsenddispls,mpi_double_precision,&
4124 segments_recv_round,roundrecvcounts,roundrecvdispls,&
4127 if (totalroundseg>0)
then
4128 segments_recv(:,recvfill+1:recvfill+totalroundseg)=segments_recv_round(:,1:totalroundseg)
4129 recvfill=recvfill+totalroundseg
4131 deallocate(segments_recv_round)
4134 if (recvfill/=totalseg)
call mpistop(
"cart_dda thick segmented receive mismatch")
4136 if (totalseg>0)
then
4137 allocate(idx(totalseg))
4138 bucketcounts(1:npixbatch)=0
4141 ipix=nint(segments_recv(1,is))
4143 ilocal=ipix-ipixstart+1
4144 bucketcounts(ilocal)=bucketcounts(ilocal)+1
4150 do ilocal=1,npixbatch
4151 bucketoffsets(ilocal+1)=bucketoffsets(ilocal)+bucketcounts(ilocal)
4153 bucketfill(1:npixbatch)=bucketoffsets(1:npixbatch)
4156 ipix=nint(segments_recv(1,is))
4158 ilocal=ipix-ipixstart+1
4159 idx(bucketfill(ilocal))=is
4160 bucketfill(ilocal)=bucketfill(ilocal)+1
4165 do ipix=ipixstart,ipixend
4167 ilocal=ipix-ipixstart+1
4168 nidx=bucketcounts(ilocal)
4170 profile_local(5)=profile_local(5)+dble(nidx)*dble(nidx)
4172 ixglobal=1+mod(ipix-1,nxif1)
4173 iyglobal=1+(ipix-1)/nxif1
4174 do iseg=bucketoffsets(ilocal),bucketoffsets(ilocal+1)-1
4176 euvthin(ixglobal,iyglobal)=euvthin(ixglobal,iyglobal)+segments_recv(3,is)
4178 euv(ixglobal,iyglobal)=euv(ixglobal,iyglobal)+atten*segments_recv(3,is)
4179 dpl(ixglobal,iyglobal)=dpl(ixglobal,iyglobal)+atten*segments_recv(5,is)
4180 tau(ixglobal,iyglobal)=tau(ixglobal,iyglobal)+max(zero,segments_recv(4,is))
4187 deallocate(segments_send,segments_recv)
4188 if (
allocated(segments))
deallocate(segments)
4192 do iigrid=1,igridstail
4193 if (
allocated(cache(iigrid)%source))
deallocate(cache(iigrid)%source)
4194 if (
allocated(cache(iigrid)%opacity))
deallocate(cache(iigrid)%opacity)
4195 if (
allocated(cache(iigrid)%sourcev))
deallocate(cache(iigrid)%sourcev)
4196 if (
allocated(cache(iigrid)%xface1))
deallocate(cache(iigrid)%xface1)
4197 if (
allocated(cache(iigrid)%xface2))
deallocate(cache(iigrid)%xface2)
4198 if (
allocated(cache(iigrid)%xface3))
deallocate(cache(iigrid)%xface3)
4201 deallocate(sendcounts,recvcounts,senddispls,recvdispls,roundsendcounts,roundrecvcounts,&
4202 roundsenddispls,roundrecvdispls,ownersegcounts,owneroffsets,bucketcounts,&
4203 bucketoffsets,bucketfill)
4204 allocate(image_reduce(nxif1,nxif2))
4205 call mpi_allreduce(euv,image_reduce,nxif1*nxif2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4207 call mpi_allreduce(dpl,image_reduce,nxif1*nxif2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4209 call mpi_allreduce(tau,image_reduce,nxif1*nxif2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4211 call mpi_allreduce(euvthin,image_reduce,nxif1*nxif2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4212 euvthin=image_reduce
4213 deallocate(image_reduce)
4214 call mpi_allreduce(profile_local,profile_global,5,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4216 write(*,
'(a,5(es12.5,1x))') &
4217 ' cart_dda thick profile: ',profile_global
4224 integer,
intent(in) :: nXIF1,nXIF2
4226 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
4227 double precision,
intent(out) :: Tau(nXIF1,nXIF2),EUVthin(nXIF1,nXIF2)
4229 integer :: ixO^L,ixO^D,ixI^L,ix^D
4230 integer :: iigrid,igrid,levelg,rft,direction_LOS,nLOS,numSeg,nLayerVars,nLayerSeg
4231 integer :: ixP1,ixP2,ixL,iSub1,iSub2,relL
4232 integer :: nLosBatch,nBatch,iBatch,ixLstart,ixLend,ixLgridStart,ixLgridEnd
4233 integer :: ixPmin1,ixPmin2
4234 double precision :: ds_cm,jds,kds,jvds,atten,layerBytes,targetBytes
4235 double precision,
allocatable :: rho(:^D&)
4236 double precision,
allocatable :: layer_ds(:,:,:,:),layer_all(:,:,:,:)
4251 if (nxif1>huge(numseg)/max(1,nxif2) .or. nxif1*nxif2>huge(numseg)/nlayervars)
then
4252 call mpistop(
"thick EUV layer buffer is too large for one MPI reduction")
4254 nlayerseg=nxif1*nxif2*nlayervars
4255 targetbytes=256.d0*1024.d0*1024.d0
4256 layerbytes=dble(nlayerseg)*8.d0*2.d0
4257 nlosbatch=max(1,min(16,int(targetbytes/max(one,layerbytes))))
4258 if (nlayerseg>huge(numseg)/nlosbatch)
then
4259 call mpistop(
"thick EUV batched layer buffer is too large for one MPI reduction")
4262 allocate(cache(igridstail))
4263 do iigrid=1,igridstail; igrid=igrids(iigrid);
4264 ^d&ixomin^d=ixmlo^d\
4265 ^d&ixomax^d=ixmhi^d\
4266 ^d&iximin^d=
ixglo^d\
4267 ^d&iximax^d=
ixghi^d\
4269 cache(iigrid)%igrid=igrid
4270 levelg=ps(igrid)%level
4272 cache(iigrid)%level=levelg
4273 cache(iigrid)%rft=rft
4275 select case(direction_los)
4277 cache(iigrid)%los_min=(
node(pig1_,igrid)-1)*rft*block_nx1+1
4278 cache(iigrid)%los_max=
node(pig1_,igrid)*rft*block_nx1
4280 cache(iigrid)%los_min=(
node(pig2_,igrid)-1)*rft*block_nx2+1
4281 cache(iigrid)%los_max=
node(pig2_,igrid)*rft*block_nx2
4283 cache(iigrid)%los_min=(
node(pig3_,igrid)-1)*rft*block_nx3+1
4284 cache(iigrid)%los_max=
node(pig3_,igrid)*rft*block_nx3
4287 allocate(cache(iigrid)%source(ixi^s),cache(iigrid)%opacity(ixi^s),&
4288 cache(iigrid)%sourcev(ixi^s),rho(ixi^s))
4289 cache(iigrid)%source=zero
4290 cache(iigrid)%opacity=zero
4291 cache(iigrid)%sourcev=zero
4294 cache(iigrid)%source,cache(iigrid)%opacity)
4296 call get_euv(
wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,cache(iigrid)%source)
4299 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
4300 cache(iigrid)%sourcev(ixo^s)=cache(iigrid)%source(ixo^s)*&
4301 (-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s))
4307 allocate(layer_ds(nxif1,nxif2,nlayervars,nlosbatch),layer_all(nxif1,nxif2,nlayervars,nlosbatch))
4313 do ixlstart=1,nlos,nlosbatch
4314 ixlend=min(nlos,ixlstart+nlosbatch-1)
4315 nbatch=ixlend-ixlstart+1
4316 layer_ds(:,:,:,1:nbatch)=zero
4318 do iigrid=1,igridstail
4319 ixlgridstart=max(ixlstart,cache(iigrid)%los_min)
4320 ixlgridend=min(ixlend,cache(iigrid)%los_max)
4321 if (ixlgridstart>ixlgridend) cycle
4322 igrid=cache(iigrid)%igrid
4323 rft=cache(iigrid)%rft
4324 ^d&ixomin^d=ixmlo^d\
4325 ^d&ixomax^d=ixmhi^d\
4326 ^d&iximin^d=
ixglo^d\
4327 ^d&iximax^d=
ixghi^d\
4329 do ixl=ixlgridstart,ixlgridend
4330 ibatch=ixl-ixlstart+1
4331 rell=ixl-cache(iigrid)%los_min
4333 select case(direction_los)
4335 ix1=ixomin1+rell/rft
4336 do ix2=ixomin2,ixomax2
4337 ixpmin1=(
node(pig2_,igrid)-1)*rft*block_nx2+(ix2-ixomin2)*rft+1
4338 do ix3=ixomin3,ixomax3
4339 ixpmin2=(
node(pig3_,igrid)-1)*rft*block_nx3+(ix3-ixomin3)*rft+1
4342 jds=cache(iigrid)%source(ix^d)*ds_cm
4343 kds=cache(iigrid)%opacity(ix^d)*ds_cm
4344 jvds=cache(iigrid)%sourcev(ix^d)*ds_cm
4349 layer_ds(ixp1,ixp2,1,ibatch)=layer_ds(ixp1,ixp2,1,ibatch)+jds
4350 layer_ds(ixp1,ixp2,2,ibatch)=layer_ds(ixp1,ixp2,2,ibatch)+kds
4351 layer_ds(ixp1,ixp2,3,ibatch)=layer_ds(ixp1,ixp2,3,ibatch)+jvds
4357 ix2=ixomin2+rell/rft
4358 do ix3=ixomin3,ixomax3
4359 ixpmin1=(
node(pig3_,igrid)-1)*rft*block_nx3+(ix3-ixomin3)*rft+1
4360 do ix1=ixomin1,ixomax1
4361 ixpmin2=(
node(pig1_,igrid)-1)*rft*block_nx1+(ix1-ixomin1)*rft+1
4364 jds=cache(iigrid)%source(ix^d)*ds_cm
4365 kds=cache(iigrid)%opacity(ix^d)*ds_cm
4366 jvds=cache(iigrid)%sourcev(ix^d)*ds_cm
4371 layer_ds(ixp1,ixp2,1,ibatch)=layer_ds(ixp1,ixp2,1,ibatch)+jds
4372 layer_ds(ixp1,ixp2,2,ibatch)=layer_ds(ixp1,ixp2,2,ibatch)+kds
4373 layer_ds(ixp1,ixp2,3,ibatch)=layer_ds(ixp1,ixp2,3,ibatch)+jvds
4379 ix3=ixomin3+rell/rft
4380 do ix1=ixomin1,ixomax1
4381 ixpmin1=(
node(pig1_,igrid)-1)*rft*block_nx1+(ix1-ixomin1)*rft+1
4382 do ix2=ixomin2,ixomax2
4383 ixpmin2=(
node(pig2_,igrid)-1)*rft*block_nx2+(ix2-ixomin2)*rft+1
4386 jds=cache(iigrid)%source(ix^d)*ds_cm
4387 kds=cache(iigrid)%opacity(ix^d)*ds_cm
4388 jvds=cache(iigrid)%sourcev(ix^d)*ds_cm
4393 layer_ds(ixp1,ixp2,1,ibatch)=layer_ds(ixp1,ixp2,1,ibatch)+jds
4394 layer_ds(ixp1,ixp2,2,ibatch)=layer_ds(ixp1,ixp2,2,ibatch)+kds
4395 layer_ds(ixp1,ixp2,3,ibatch)=layer_ds(ixp1,ixp2,3,ibatch)+jvds
4404 numseg=nlayerseg*nbatch
4405 call mpi_allreduce(layer_ds,layer_all,numseg,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4411 euvthin(ixp1,ixp2)=euvthin(ixp1,ixp2)+layer_all(ixp1,ixp2,1,ibatch)
4413 euv(ixp1,ixp2)=euv(ixp1,ixp2)+atten*layer_all(ixp1,ixp2,1,ibatch)
4414 dpl(ixp1,ixp2)=dpl(ixp1,ixp2)+atten*layer_all(ixp1,ixp2,3,ibatch)
4415 tau(ixp1,ixp2)=tau(ixp1,ixp2)+layer_all(ixp1,ixp2,2,ibatch)
4422 do iigrid=1,igridstail
4423 if (
allocated(cache(iigrid)%source))
deallocate(cache(iigrid)%source)
4424 if (
allocated(cache(iigrid)%opacity))
deallocate(cache(iigrid)%opacity)
4425 if (
allocated(cache(iigrid)%sourcev))
deallocate(cache(iigrid)%sourcev)
4427 deallocate(cache,layer_ds,layer_all)
4436 integer,
intent(in) :: ixI^L, ixO^L
4437 double precision,
intent(in) :: x(ixI^S,1:ndim),dx(ixI^S,1:ndim)
4438 double precision,
allocatable,
intent(out) :: rface(:),thetaface(:),phiface(:)
4439 integer :: ix1,ix2,ix3
4441 allocate(rface(ixomin1:ixomax1+1),thetaface(ixomin2:ixomax2+1),phiface(ixomin3:ixomax3+1))
4442 do ix1=ixomin1,ixomax1
4443 rface(ix1)=x(ix1,ixomin2,ixomin3,1)-half*dx(ix1,ixomin2,ixomin3,1)
4445 rface(ixomax1+1)=x(ixomax1,ixomin2,ixomin3,1)+half*dx(ixomax1,ixomin2,ixomin3,1)
4446 do ix2=ixomin2,ixomax2
4447 thetaface(ix2)=x(ixomin1,ix2,ixomin3,2)-half*dx(ixomin1,ix2,ixomin3,2)
4449 thetaface(ixomax2+1)=x(ixomin1,ixomax2,ixomin3,2)+half*dx(ixomin1,ixomax2,ixomin3,2)
4450 do ix3=ixomin3,ixomax3
4451 phiface(ix3)=x(ixomin1,ixomin2,ix3,3)-half*dx(ixomin1,ixomin2,ix3,3)
4453 phiface(ixomax3+1)=x(ixomin1,ixomin2,ixomax3,3)+half*dx(ixomin1,ixomin2,ixomax3,3)
4457 double precision,
allocatable,
intent(inout) :: tvals(:)
4458 integer,
intent(inout) :: nt,capacity
4459 double precision,
intent(in) :: t
4461 double precision,
allocatable :: tmp(:)
4463 if (t /= t .or. abs(t)>1.d90)
return
4464 if (.not.
allocated(tvals))
then
4466 allocate(tvals(capacity))
4467 else if (nt>=capacity)
then
4468 allocate(tmp(capacity))
4471 allocate(tvals(2*capacity))
4472 tvals(1:capacity)=tmp
4481 double precision,
intent(inout) :: tvals(:)
4482 integer,
intent(inout) :: nt
4483 double precision,
intent(in) :: t
4485 if (t /= t .or. abs(t)>1.d90)
return
4486 if (nt>=
size(tvals))
return
4492 double precision,
intent(inout) :: tvals(:)
4493 integer,
intent(inout) :: nt
4496 double precision :: key,epsT
4502 do while (j>=1 .and. tvals(j)>key)
4508 epst=max(1.d-12,1.d-10*max(one,abs(tvals(nt)-tvals(1))))
4511 if (abs(tvals(i)-tvals(nout))>epst)
then
4513 tvals(nout)=tvals(i)
4520 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),rface
4521 double precision,
allocatable,
intent(inout) :: tvals(:)
4522 integer,
intent(inout) :: nt,capacity
4524 double precision :: aa,bb,cc,disc,root
4527 bb=2.d0*sum(ray_origin*ray_dir)
4528 cc=sum(ray_origin**2)-rface**2
4529 disc=bb**2-4.d0*aa*cc
4530 if (disc<zero)
return
4531 root=sqrt(max(zero,disc))
4532 call sph_add_t(tvals,nt,capacity,(-bb-root)/(2.d0*aa))
4533 call sph_add_t(tvals,nt,capacity,(-bb+root)/(2.d0*aa))
4537 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),thetaface
4538 double precision,
allocatable,
intent(inout) :: tvals(:)
4539 integer,
intent(inout) :: nt,capacity
4541 double precision :: cth,aa,bb,cc,disc,root
4547 if (abs(cth)<1.d-12)
then
4548 if (abs(ray_dir(3))>1.d-14) &
4549 call sph_add_t(tvals,nt,capacity,-ray_origin(3)/ray_dir(3))
4552 aa=ray_dir(3)**2-cth**2*sum(ray_dir**2)
4553 bb=2.d0*(ray_origin(3)*ray_dir(3)-cth**2*sum(ray_origin*ray_dir))
4554 cc=ray_origin(3)**2-cth**2*sum(ray_origin**2)
4555 if (abs(aa)<1.d-14)
then
4556 if (abs(bb)>1.d-14)
call sph_add_t(tvals,nt,capacity,-cc/bb)
4559 disc=bb**2-4.d0*aa*cc
4560 if (disc<zero)
return
4561 root=sqrt(max(zero,disc))
4562 call sph_add_t(tvals,nt,capacity,(-bb-root)/(2.d0*aa))
4563 call sph_add_t(tvals,nt,capacity,(-bb+root)/(2.d0*aa))
4567 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),phiface
4568 double precision,
allocatable,
intent(inout) :: tvals(:)
4569 integer,
intent(inout) :: nt,capacity
4571 double precision :: normal(1:3),denom,numer
4573 normal(1)=-sin(phiface)
4574 normal(2)=cos(phiface)
4576 denom=sum(normal*ray_dir)
4577 if (abs(denom)<1.d-14)
return
4578 numer=sum(normal*ray_origin)
4579 call sph_add_t(tvals,nt,capacity,-numer/denom)
4583 double precision,
intent(in) :: pos(1:3)
4584 double precision,
intent(out) :: sph(1:3)
4586 sph(1)=sqrt(sum(pos**2))
4587 if (sph(1)>zero)
then
4588 sph(2)=acos(max(-one,min(one,pos(3)/sph(1))))
4592 sph(3)=atan2(pos(2),pos(1))
4596 integer,
intent(in) :: imin,imax
4597 double precision,
intent(in) ::
value,faces(imin:imax+1)
4599 integer :: ilo,ihi,imid
4602 if (
value<faces(imin)-1.d-12 .or.
value>faces(imax+1)+1.d-12)
return
4603 if (
value<=faces(imin))
then
4607 if (
value>=faces(imax+1))
then
4614 do while (ihi-ilo>1)
4616 if (
value>=faces(imid))
then
4622 idx=min(imax,max(imin,ilo))
4626 integer,
intent(in) :: imin,imax
4627 double precision,
intent(in) ::
value,faces(imin:imax+1)
4629 integer :: ilo,ihi,imid
4632 if (
value>faces(imin)+1.d-12 .or.
value<faces(imax+1)-1.d-12)
return
4633 if (
value>=faces(imin))
then
4637 if (
value<=faces(imax+1))
then
4644 do while (ihi-ilo>1)
4646 if (
value<=faces(imid))
then
4652 idx=min(imax,max(imin,ilo))
4656 double precision,
intent(in) :: pos(1:3)
4657 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4658 phiface(ixOmin3:ixOmax3+1)
4659 integer,
intent(in) :: ixO^L
4660 integer,
intent(out) :: ix1,ix2,ix3
4661 logical,
intent(out) :: inside
4663 double precision :: sph(1:3),phi
4667 if (phi<phiface(ixomin3)-1.d-12) phi=phi+2.d0*dpi
4668 if (phi>phiface(ixomax3+1)+1.d-12) phi=phi-2.d0*dpi
4672 inside=ix1>0 .and. ix2>0 .and. ix3>0
4677 double precision,
intent(in) :: pos(1:3),ximg1,ximg2
4679 double precision :: dotp,rc,rthick,rloc
4682 rc=sqrt(ximg1**2+ximg2**2)
4683 rloc=sqrt(sum(pos**2))
4686 if (dotp>=
zero)
then
4687 if (rc<=rthick) visible=.false.
4689 if (rloc<=rthick) visible=.false.
4694 ixPmin1,ixPmax1,ixPmin2,ixPmax2,has_pixels)
4695 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4696 phiface(ixOmin3:ixOmax3+1)
4697 integer,
intent(in) :: ixO^L,nXI1,nXI2
4698 double precision,
intent(in) :: xI1(nXI1),xI2(nXI2),dxI
4699 integer,
intent(out) :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
4700 logical,
intent(out) :: has_pixels
4702 integer,
parameter :: nsample=5
4704 double precision :: sph(1:3),xcent(1:2)
4705 double precision :: xmin1,xmax1,xmin2,xmax2
4706 double precision :: wr,wt,wp,pad
4714 wr=dble(ir)/dble(nsample-1)
4715 sph(1)=(one-wr)*rface(ixomin1)+wr*rface(ixomax1+1)
4717 wt=dble(it)/dble(nsample-1)
4718 sph(2)=(one-wt)*thetaface(ixomin2)+wt*thetaface(ixomax2+1)
4720 wp=dble(ip)/dble(nsample-1)
4721 if (ir/=0 .and. ir/=nsample-1 .and. it/=0 .and. it/=nsample-1 .and. &
4722 ip/=0 .and. ip/=nsample-1) cycle
4723 sph(3)=(one-wp)*phiface(ixomin3)+wp*phiface(ixomax3+1)
4725 xmin1=min(xmin1,xcent(1))
4726 xmax1=max(xmax1,xcent(1))
4727 xmin2=min(xmin2,xcent(2))
4728 xmax2=max(xmax2,xcent(2))
4737 ixpmin1=max(1,floor((xmin1-(xi1(1)-half*dxi))/dxi)+1)
4738 ixpmax1=min(nxi1,ceiling((xmax1-(xi1(1)-half*dxi))/dxi))
4739 ixpmin2=max(1,floor((xmin2-(xi2(1)-half*dxi))/dxi)+1)
4740 ixpmax2=min(nxi2,ceiling((xmax2-(xi2(1)-half*dxi))/dxi))
4741 has_pixels=ixpmin1<=ixpmax1 .and. ixpmin2<=ixpmax2
4745 rface,thetaface,phiface,EUVp)
4746 integer,
intent(in) :: ixI^L,ixO^L
4747 double precision,
intent(in) :: source(ixI^S),ray_origin(1:3),ximg1,ximg2
4748 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4749 phiface(ixOmin3:ixOmax3+1)
4750 double precision,
intent(inout) :: EUVp
4752 integer :: nt,capacity,i,ix^D
4753 double precision,
allocatable :: tvals(:)
4754 double precision :: posMid(1:3),ds_cm,tMid,t0,t1
4759 do ix1=ixomin1,ixomax1+1
4762 do ix2=ixomin2,ixomax2+1
4765 do ix3=ixomin3,ixomax3+1
4769 if (
allocated(tvals))
deallocate(tvals)
4779 posmid=ray_origin+tmid*
vec_los
4781 call sph_locate_cell(posmid,rface,thetaface,phiface,ixo^l,ix1,ix2,ix3,inside)
4782 if (.not. inside) cycle
4783 ds_cm=(t1-t0)*unit_length
4784 if (si_unit) ds_cm=ds_cm*1.d2
4785 euvp=euvp+source(ix^d)*ds_cm
4791 pixel_id,ray_origin,ximg1,ximg2,&
4792 rface,thetaface,phiface,rface2,&
4793 theta_cos,phi_sin,phi_cos,&
4794 segments,nseg,capacity)
4797 integer,
intent(in) :: ixI^L,ixO^L,pixel_id
4798 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S)
4799 double precision,
intent(in) :: ray_origin(1:3),ximg1,ximg2
4800 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4801 phiface(ixOmin3:ixOmax3+1)
4802 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4803 phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
4804 double precision,
allocatable,
intent(inout) :: segments(:,:)
4805 integer,
intent(inout) :: nseg,capacity
4807 integer :: nt,i,ix^D
4808 double precision :: tvals(2*(ixOmax1-ixOmin1+2)+2*(ixOmax2-ixOmin2+2)+&
4809 (ixOmax3-ixOmin3+2))
4810 double precision :: posMid(1:3),ds_cm,tMid,t0,t1,jds,kds
4811 double precision :: dir2,origin2,odotd,aa,bb,cc,disc,root,cth2
4812 double precision :: denom,numer,r2,mu,phi,dotp,rthick2,rc2
4817 origin2=sum(ray_origin**2)
4819 rthick2=(r_opt_thick*
const_rsun/unit_length)**2
4820 rc2=ximg1**2+ximg2**2
4822 do ix1=ixomin1,ixomax1+1
4825 cc=origin2-rface2(ix1)
4826 disc=bb**2-4.d0*aa*cc
4827 if (disc>=
zero)
then
4828 root=sqrt(max(
zero,disc))
4833 do ix2=ixomin2,ixomax2+1
4834 if (abs(theta_cos(ix2))<1.d-12)
then
4839 cth2=theta_cos(ix2)**2
4841 bb=2.d0*(ray_origin(3)*
vec_los(3)-cth2*odotd)
4842 cc=ray_origin(3)**2-cth2*origin2
4843 if (abs(aa)<1.d-14)
then
4846 disc=bb**2-4.d0*aa*cc
4847 if (disc>=
zero)
then
4848 root=sqrt(max(
zero,disc))
4854 do ix3=ixomin3,ixomax3+1
4856 if (abs(denom)>=1.d-14)
then
4857 numer=-phi_sin(ix3)*ray_origin(1)+phi_cos(ix3)*ray_origin(2)
4869 posmid=ray_origin+tmid*
vec_los
4872 if (dotp>=
zero)
then
4873 if (rc2<=rthick2) cycle
4875 if (r2<=rthick2) cycle
4880 mu=posmid(3)/sqrt(r2)
4885 phi=atan2(posmid(2),posmid(1))
4886 if (phi<phiface(ixomin3)-1.d-12) phi=phi+2.d0*
dpi
4887 if (phi>phiface(ixomax3+1)+1.d-12) phi=phi-2.d0*
dpi
4889 inside=ix1>0 .and. ix2>0 .and. ix3>0
4890 if (.not. inside) cycle
4891 ds_cm=(t1-t0)*unit_length
4892 if (si_unit) ds_cm=ds_cm*1.d2
4893 jds=max(
zero,source(ix^d))*ds_cm
4894 kds=max(
zero,opacity(ix^d))*ds_cm
4900 double precision,
intent(in) :: pos(1:3)
4901 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4902 phiface(ixOmin3:ixOmax3+1)
4903 integer,
intent(in) :: ixO^L
4904 integer,
intent(out) :: ix1,ix2,ix3
4905 logical,
intent(out) :: inside
4907 double precision :: r2,mu,phi
4917 phi=atan2(pos(2),pos(1))
4918 if (phi<phiface(ixomin3)-1.d-12) phi=phi+2.d0*dpi
4919 if (phi>phiface(ixomax3+1)+1.d-12) phi=phi-2.d0*dpi
4921 inside=ix1>0 .and. ix2>0 .and. ix3>0
4925 double precision,
intent(in) :: t,tNow,tExit,epsRay
4926 double precision,
intent(inout) :: tNext
4927 logical,
intent(inout) :: found
4929 if (t>tnow+epsray .and. t<=texit+epsray .and. t<tnext)
then
4936 double precision,
intent(in) :: t,theta_face_cos,ray_origin(1:3),tNow,tExit,epsRay
4937 double precision,
intent(inout) :: tNext
4938 logical,
intent(inout) :: found
4940 double precision :: pos(1:3),r2
4942 if (t<=tnow+epsray .or. t>texit+epsray .or. t>=tnext)
return
4945 if (r2<=zero)
return
4946 if (theta_face_cos>1.d-12 .and. pos(3)<-1.d-10)
return
4947 if (theta_face_cos<-1.d-12 .and. pos(3)>1.d-10)
return
4948 if (abs(pos(3)**2-theta_face_cos**2*r2)>1.d-6*max(one,r2))
return
4954 double precision,
intent(in) :: t,phi_face_sin,phi_face_cos,ray_origin(1:3),tNow,tExit,epsRay
4955 double precision,
intent(inout) :: tNext
4956 logical,
intent(inout) :: found
4958 double precision :: pos(1:3),radialDot
4960 if (t<=tnow+epsray .or. t>texit+epsray .or. t>=tnext)
return
4962 radialdot=phi_face_cos*pos(1)+phi_face_sin*pos(2)
4963 if (radialdot<-1.d-10)
return
4969 ix1,ix2,ix3,tNow,tExit,epsRay,tNext,found)
4970 double precision,
intent(in) :: ray_origin(1:3)
4971 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4972 phiface(ixOmin3:ixOmax3+1)
4973 double precision,
intent(in) :: phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
4974 integer,
intent(in) :: ixO^L,ix1,ix2,ix3
4975 double precision,
intent(in) :: tNow,tExit,epsRay
4976 double precision,
intent(out) :: tNext
4977 logical,
intent(out) :: found
4980 double precision :: dir2,origin2,odotd,aa,bb,cc,disc,root,cth2,denom,numer
4985 origin2=sum(ray_origin**2)
4991 cc=origin2-rface2(iface)
4992 disc=bb**2-4.d0*aa*cc
4993 if (disc>=zero)
then
4994 root=sqrt(max(zero,disc))
5001 if (abs(theta_cos(iface))<1.d-12)
then
5003 -ray_origin(3)/
vec_los(3),theta_cos(iface),ray_origin,tnow,texit,epsray,&
5007 cth2=theta_cos(iface)**2
5009 bb=2.d0*(ray_origin(3)*
vec_los(3)-cth2*odotd)
5010 cc=ray_origin(3)**2-cth2*origin2
5011 if (abs(aa)<1.d-14)
then
5013 ray_origin,tnow,texit,epsray,tnext,found)
5015 disc=bb**2-4.d0*aa*cc
5016 if (disc>=zero)
then
5017 root=sqrt(max(zero,disc))
5019 ray_origin,tnow,texit,epsray,tnext,found)
5021 ray_origin,tnow,texit,epsray,tnext,found)
5028 if (abs(denom)>=1.d-14)
then
5029 numer=-phi_sin(iface)*ray_origin(1)+phi_cos(iface)*ray_origin(2)
5031 tnow,texit,epsray,tnext,found)
5037 ray_origin,ximg1,ximg2,rface2,theta_cos,phiface,&
5039 t_enter,t_exit,segments,nseg,capacity,ok)
5042 integer,
intent(in) :: ixI^L,ixO^L,pixel_id
5043 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S)
5044 double precision,
intent(in) :: ray_origin(1:3),ximg1,ximg2
5045 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
5046 phiface(ixOmin3:ixOmax3+1)
5047 double precision,
intent(in) :: phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
5048 double precision,
intent(in) :: t_enter,t_exit
5049 double precision,
allocatable,
intent(inout) :: segments(:,:)
5050 integer,
intent(inout) :: nseg,capacity
5051 logical,
intent(out) :: ok
5053 integer :: ix^D,nstep,maxSteps
5054 double precision :: tNow,tNext,tEnd,tMid,epsRay,ds_cm,jds,kds
5055 double precision :: pos(1:3),r2,dotp,rthick2,rc2
5056 logical :: inside,found
5059 if (t_exit<=t_enter)
return
5060 epsray=max(1.d-12,1.d-10*max(
one,abs(t_exit-t_enter)))
5061 pos=ray_origin+(t_enter+epsray)*
vec_los
5063 if (.not. inside)
then
5068 rthick2=(r_opt_thick*
const_rsun/unit_length)**2
5069 rc2=ximg1**2+ximg2**2
5072 maxsteps=8*((ixomax1-ixomin1+1)+(ixomax2-ixomin2+1)+(ixomax3-ixomin3+1)+3)
5074 do while (tnow<t_exit-epsray)
5076 ix1,ix2,ix3,tnow,t_exit,epsray,tnext,found)
5077 if (.not. found)
then
5081 tend=min(tnext,t_exit)
5083 tmid=
half*(tnow+tend)
5087 if (.not. ((dotp>=
zero .and. rc2<=rthick2) .or. &
5088 (dotp<
zero .and. r2<=rthick2)))
then
5089 ds_cm=(tend-tnow)*unit_length
5090 if (si_unit) ds_cm=ds_cm*1.d2
5091 jds=max(
zero,source(ix^d))*ds_cm
5092 kds=max(
zero,opacity(ix^d))*ds_cm
5097 if (tnow>=t_exit-epsray)
exit
5098 pos=ray_origin+(tnow+epsray)*
vec_los
5100 if (.not. inside)
then
5105 if (nstep>maxsteps)
then
5113 pixel_id,ray_origin,ximg1,ximg2,&
5114 rface,thetaface,phiface,rface2,&
5115 theta_cos,phi_sin,phi_cos,&
5116 segments,nseg,capacity,fallback)
5117 integer,
intent(in) :: ixI^L,ixO^L,pixel_id
5118 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S)
5119 double precision,
intent(in) :: ray_origin(1:3),ximg1,ximg2
5120 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
5121 phiface(ixOmin3:ixOmax3+1)
5122 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
5123 phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
5124 double precision,
allocatable,
intent(inout) :: segments(:,:)
5125 integer,
intent(inout) :: nseg,capacity
5126 logical,
intent(out) :: fallback
5128 integer :: nt,i,nsegStart,ix^D
5129 double precision :: tvals(12),posMid(1:3),t0,t1,tMid
5130 double precision :: dir2,origin2,odotd,aa,bb,cc,disc,root,cth2,denom,numer
5131 logical :: inside,ok
5137 origin2=sum(ray_origin**2)
5140 do ix1=ixomin1,ixomax1+1,ixomax1-ixomin1+1
5143 cc=origin2-rface2(ix1)
5144 disc=bb**2-4.d0*aa*cc
5145 if (disc>=zero)
then
5146 root=sqrt(max(zero,disc))
5151 do ix2=ixomin2,ixomax2+1,ixomax2-ixomin2+1
5152 if (abs(theta_cos(ix2))<1.d-12)
then
5157 cth2=theta_cos(ix2)**2
5159 bb=2.d0*(ray_origin(3)*
vec_los(3)-cth2*odotd)
5160 cc=ray_origin(3)**2-cth2*origin2
5161 if (abs(aa)<1.d-14)
then
5164 disc=bb**2-4.d0*aa*cc
5165 if (disc>=zero)
then
5166 root=sqrt(max(zero,disc))
5172 do ix3=ixomin3,ixomax3+1,ixomax3-ixomin3+1
5174 if (abs(denom)>=1.d-14)
then
5175 numer=-phi_sin(ix3)*ray_origin(1)+phi_cos(ix3)*ray_origin(2)
5188 posmid=ray_origin+tmid*
vec_los
5190 if (.not. inside) cycle
5192 ray_origin,ximg1,ximg2,rface2,theta_cos,phiface,phi_sin,phi_cos,&
5193 t0,t1,segments,nseg,capacity,ok)
5198 pixel_id,ray_origin,ximg1,ximg2,rface,thetaface,phiface,rface2,&
5199 theta_cos,phi_sin,phi_cos,segments,nseg,capacity)
5209 integer,
intent(in) :: numXI1,numXI2
5210 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2),dxI
5212 double precision,
intent(inout) :: EM(numXI1,numXI2)
5214 integer :: ixO^L,ixI^L,ix^D
5215 integer :: iigrid,igrid,ixP1,ixP2,ixPmin1,ixPmax1,ixPmin2,ixPmax2
5216 integer :: iseg,nseg,capacity,sphDdaFallbackLocal,sphDdaFallbackGlobal
5217 double precision :: ray_origin(1:3),profile_local(3),profile_global(3)
5218 double precision,
allocatable :: source(:^D&)
5219 double precision,
allocatable :: rface(:),thetaface(:),phiface(:)
5220 double precision,
allocatable :: rface2(:),theta_cos(:),phi_sin(:),phi_cos(:)
5221 double precision,
allocatable :: segments(:,:)
5222 logical :: has_pixels,ddaFallback
5225 sphddafallbacklocal=0
5226 do iigrid=1,igridstail; igrid=igrids(iigrid);
5227 ^d&ixomin^d=ixmlo^d\
5228 ^d&ixomax^d=ixmhi^d\
5229 ^d&iximin^d=
ixglo^d\
5230 ^d&iximax^d=
ixghi^d\
5232 allocate(source(ixi^s))
5237 allocate(rface2(ixomin1:ixomax1+1),theta_cos(ixomin2:ixomax2+1),&
5238 phi_sin(ixomin3:ixomax3+1),phi_cos(ixomin3:ixomax3+1))
5240 theta_cos=cos(thetaface)
5241 phi_sin=sin(phiface)
5242 phi_cos=cos(phiface)
5245 ixpmin1,ixpmax1,ixpmin2,ixpmax2,has_pixels)
5246 if (has_pixels)
then
5248 do ixp1=ixpmin1,ixpmax1
5249 do ixp2=ixpmin2,ixpmax2
5251 profile_local(1)=profile_local(1)+one
5255 1,ray_origin,xi1(ixp1),xi2(ixp2),rface,thetaface,phiface,&
5256 rface2,theta_cos,phi_sin,phi_cos,segments,nseg,capacity,ddafallback)
5257 if (ddafallback) sphddafallbacklocal=sphddafallbacklocal+1
5259 em(ixp1,ixp2)=em(ixp1,ixp2)+segments(3,iseg)
5263 rface,thetaface,phiface,em(ixp1,ixp2))
5267 profile_local(2)=profile_local(2)+dble((ixpmax1-ixpmin1+1)*(ixpmax2-ixpmin2+1))
5269 if (
allocated(segments))
deallocate(segments)
5270 profile_local(3)=profile_local(3)+one
5271 if (
allocated(rface2))
deallocate(rface2)
5272 if (
allocated(theta_cos))
deallocate(theta_cos)
5273 if (
allocated(phi_sin))
deallocate(phi_sin)
5274 if (
allocated(phi_cos))
deallocate(phi_cos)
5275 deallocate(source,rface,thetaface,phiface)
5277 call mpi_allreduce(profile_local,profile_global,3,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5278 call mpi_allreduce(sphddafallbacklocal,sphddafallbackglobal,1,mpi_integer,mpi_sum,
icomm,
ierrmpi)
5281 write(*,
'(a,3(es12.5,1x))')
' sph_dda thin profile rays pixels blocks: ',profile_global
5282 write(*,
'(a,i0)')
' sph_dda thin fallback rays: ',sphddafallbackglobal
5284 write(*,
'(a,3(es12.5,1x))')
' sph_intersection thin profile rays pixels blocks: ',profile_global
5292 integer,
intent(in) :: numXI1,numXI2
5293 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2),dxI
5295 double precision,
intent(out) :: EUV(numXI1,numXI2),Tau(numXI1,numXI2),EUVthin(numXI1,numXI2)
5297 integer,
parameter :: nSegVars=5
5298 integer :: ixO^L,ixI^L,ix^D
5299 integer :: iigrid,igrid,ixP1,ixP2,ipix,ipixStart,ipixEnd,nPixBatch,pixel_id
5300 integer :: nseg,capacity,totalCount,totalSeg,ipe,is,iseg,nidx,owner,isegDest,nsegBefore
5301 integer :: ixGlobal,iyGlobal,ixPmin1,ixPmax1,ixPmin2,ixPmax2,iFirst,iLast,iLocal
5302 integer :: nPixBatchTarget
5303 integer :: maxSegBatchTarget,maxSegCommTarget,maxNsegBatch,nPixTotal
5304 integer :: maxOwnerSegCount,maxOwnerSegCountLocal,segOffset,recvFill,totalRoundCount,totalRoundSeg
5305 integer :: sphDdaFallbackLocal,sphDdaFallbackGlobal
5306 integer,
allocatable :: sendCounts(:),recvCounts(:),sendDispls(:),recvDispls(:)
5307 integer,
allocatable :: roundSendCounts(:),roundRecvCounts(:)
5308 integer,
allocatable :: roundSendDispls(:),roundRecvDispls(:)
5309 integer,
allocatable :: ownerSegCounts(:),ownerOffsets(:),idx(:)
5310 integer,
allocatable :: bucketCounts(:),bucketOffsets(:),bucketFill(:)
5311 double precision :: ray_origin(1:3),atten
5312 double precision :: profile_local(5),profile_global(5),profile_batch(5)
5313 double precision :: phys_max_local(2),phys_max_global(2)
5314 double precision :: phys_sum_local(2),phys_sum_global(2),phys_sum_batch(2)
5315 double precision,
allocatable :: segments(:,:),segments_send(:,:),segments_recv(:,:)
5316 double precision,
allocatable :: segments_recv_round(:,:)
5317 double precision,
allocatable :: image_reduce(:,:)
5318 logical :: has_pixels,batchAccepted,batchReduced,ddaFallback
5327 sphddafallbacklocal=0
5328 allocate(sendcounts(0:
npe-1),recvcounts(0:
npe-1),senddispls(0:
npe-1),recvdispls(0:
npe-1))
5329 allocate(roundsendcounts(0:
npe-1),roundrecvcounts(0:
npe-1))
5330 allocate(roundsenddispls(0:
npe-1),roundrecvdispls(0:
npe-1))
5331 allocate(ownersegcounts(0:
npe-1),owneroffsets(0:
npe-1))
5332 allocate(cache(igridstail))
5334 allocate(bucketcounts(npixbatchtarget),bucketoffsets(npixbatchtarget+1),&
5335 bucketfill(npixbatchtarget))
5337 do iigrid=1,igridstail; igrid=igrids(iigrid);
5338 ^d&ixomin^d=ixmlo^d\
5339 ^d&ixomax^d=ixmhi^d\
5340 ^d&iximin^d=
ixglo^d\
5341 ^d&iximax^d=
ixghi^d\
5343 cache(iigrid)%igrid=igrid
5344 allocate(cache(iigrid)%source(ixi^s),cache(iigrid)%opacity(ixi^s))
5345 cache(iigrid)%source=zero
5346 cache(iigrid)%opacity=zero
5347 call get_euv(
wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,cache(iigrid)%source)
5350 phys_max_local(1)=max(phys_max_local(1),maxval(cache(iigrid)%source(ixo^s)))
5351 phys_max_local(2)=max(phys_max_local(2),maxval(cache(iigrid)%opacity(ixo^s)))
5353 cache(iigrid)%rface,cache(iigrid)%thetaface,&
5354 cache(iigrid)%phiface)
5355 allocate(cache(iigrid)%rface2(ixomin1:ixomax1+1),&
5356 cache(iigrid)%theta_cos(ixomin2:ixomax2+1),&
5357 cache(iigrid)%phi_sin(ixomin3:ixomax3+1),&
5358 cache(iigrid)%phi_cos(ixomin3:ixomax3+1))
5359 cache(iigrid)%rface2=cache(iigrid)%rface**2
5360 cache(iigrid)%theta_cos=cos(cache(iigrid)%thetaface)
5361 cache(iigrid)%phi_sin=sin(cache(iigrid)%phiface)
5362 cache(iigrid)%phi_cos=cos(cache(iigrid)%phiface)
5364 cache(iigrid)%phiface,ixo^l,numxi1,numxi2,xi1,xi2,dxi,&
5365 cache(iigrid)%ixPmin1,cache(iigrid)%ixPmax1,&
5366 cache(iigrid)%ixPmin2,cache(iigrid)%ixPmax2,cache(iigrid)%has_pixels)
5369 npixtotal=numxi1*numxi2
5371 do while (ipixstart<=npixtotal)
5372 ipixend=min(numxi1*numxi2,ipixstart+npixbatchtarget-1)
5373 npixbatch=ipixend-ipixstart+1
5374 batchaccepted=.false.
5375 batchreduced=.false.
5377 do while (.not. batchaccepted)
5383 do iigrid=1,igridstail; igrid=igrids(iigrid);
5384 ^d&ixomin^d=ixmlo^d\
5385 ^d&ixomax^d=ixmhi^d\
5386 ^d&iximin^d=
ixglo^d\
5387 ^d&iximax^d=
ixghi^d\
5389 ixpmin1=cache(iigrid)%ixPmin1
5390 ixpmax1=cache(iigrid)%ixPmax1
5391 ixpmin2=cache(iigrid)%ixPmin2
5392 ixpmax2=cache(iigrid)%ixPmax2
5393 has_pixels=cache(iigrid)%has_pixels
5394 if (.not. has_pixels) cycle
5396 do ixp2=ixpmin2,ixpmax2
5397 ifirst=max(ipixstart,(ixp2-1)*numxi1+ixpmin1)
5398 ilast=min(ipixend,(ixp2-1)*numxi1+ixpmax1)
5399 if (ifirst>ilast) cycle
5400 do ipix=ifirst,ilast
5401 ixp1=1+mod(ipix-1,numxi1)
5404 profile_batch(1)=profile_batch(1)+one
5408 cache(iigrid)%opacity,pixel_id,ray_origin,xi1(ixp1),xi2(ixp2),&
5409 cache(iigrid)%rface,cache(iigrid)%thetaface,cache(iigrid)%phiface,&
5410 cache(iigrid)%rface2,cache(iigrid)%theta_cos,&
5411 cache(iigrid)%phi_sin,cache(iigrid)%phi_cos,&
5412 segments,nseg,capacity,ddafallback)
5413 if (ddafallback) sphddafallbacklocal=sphddafallbacklocal+1
5416 cache(iigrid)%opacity,pixel_id,ray_origin,xi1(ixp1),xi2(ixp2),&
5417 cache(iigrid)%rface,cache(iigrid)%thetaface,cache(iigrid)%phiface,&
5418 cache(iigrid)%rface2,cache(iigrid)%theta_cos,&
5419 cache(iigrid)%phi_sin,cache(iigrid)%phi_cos,&
5420 segments,nseg,capacity)
5422 if (nseg>nsegbefore) profile_batch(2)=profile_batch(2)+one
5423 profile_batch(3)=profile_batch(3)+dble(nseg-nsegbefore)
5424 do iseg=nsegbefore+1,nseg
5425 phys_sum_batch(1)=phys_sum_batch(1)+segments(3,iseg)
5426 phys_sum_batch(2)=phys_sum_batch(2)+segments(4,iseg)
5432 call mpi_allreduce(nseg,maxnsegbatch,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
5433 if (maxnsegbatch>maxsegbatchtarget .and. npixbatch>1)
then
5434 npixbatch=max(1,npixbatch/2)
5435 ipixend=ipixstart+npixbatch-1
5436 if (
allocated(segments))
deallocate(segments)
5439 batchaccepted=.true.
5443 profile_local=profile_local+profile_batch
5444 phys_sum_local=phys_sum_local+phys_sum_batch
5447 write(*,
'(a,3(i0,1x))')
' sph_dda thick adaptive batch: ',&
5448 ipixstart,ipixend,maxnsegbatch
5450 write(*,
'(a,3(i0,1x))')
' sph_intersection thick adaptive batch: ',&
5451 ipixstart,ipixend,maxnsegbatch
5455 if (.not.
allocated(segments))
then
5457 allocate(segments(nsegvars,capacity))
5462 ownersegcounts(owner)=ownersegcounts(owner)+1
5464 sendcounts=nsegvars*ownersegcounts
5467 senddispls(ipe)=senddispls(ipe-1)+sendcounts(ipe-1)
5470 allocate(segments_send(nsegvars,max(1,nseg)))
5474 isegdest=senddispls(owner)/nsegvars+owneroffsets(owner)+1
5475 segments_send(:,isegdest)=segments(:,is)
5476 owneroffsets(owner)=owneroffsets(owner)+1
5479 call mpi_alltoall(sendcounts,1,mpi_integer,recvcounts,1,mpi_integer,
icomm,
ierrmpi)
5482 recvdispls(ipe)=recvdispls(ipe-1)+recvcounts(ipe-1)
5484 totalcount=sum(recvcounts)
5485 totalseg=totalcount/nsegvars
5486 profile_local(4)=profile_local(4)+dble(totalcount)
5487 allocate(segments_recv(nsegvars,max(1,totalseg)))
5490 maxownersegcountlocal=maxval(ownersegcounts)
5491 call mpi_allreduce(maxownersegcountlocal,maxownersegcount,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
5492 do segoffset=0,maxownersegcount-1,maxsegcommtarget
5494 roundsenddispls=senddispls
5496 if (ownersegcounts(ipe)>segoffset)
then
5497 roundsendcounts(ipe)=nsegvars*min(maxsegcommtarget,ownersegcounts(ipe)-segoffset)
5498 roundsenddispls(ipe)=senddispls(ipe)+nsegvars*segoffset
5502 call mpi_alltoall(roundsendcounts,1,mpi_integer,roundrecvcounts,1,mpi_integer,
icomm,
ierrmpi)
5503 roundrecvdispls(0)=0
5505 roundrecvdispls(ipe)=roundrecvdispls(ipe-1)+roundrecvcounts(ipe-1)
5507 totalroundcount=sum(roundrecvcounts)
5508 totalroundseg=totalroundcount/nsegvars
5509 allocate(segments_recv_round(nsegvars,max(1,totalroundseg)))
5511 call mpi_alltoallv(segments_send,roundsendcounts,roundsenddispls,mpi_double_precision,&
5512 segments_recv_round,roundrecvcounts,roundrecvdispls,&
5515 if (totalroundseg>0)
then
5516 segments_recv(:,recvfill+1:recvfill+totalroundseg)=segments_recv_round(:,1:totalroundseg)
5517 recvfill=recvfill+totalroundseg
5519 deallocate(segments_recv_round)
5522 if (recvfill/=totalseg)
call mpistop(
"ray-segment receive mismatch")
5524 if (totalseg>0)
then
5525 allocate(idx(totalseg))
5526 bucketcounts(1:npixbatch)=0
5529 ipix=nint(segments_recv(1,is))
5531 ilocal=ipix-ipixstart+1
5532 bucketcounts(ilocal)=bucketcounts(ilocal)+1
5538 do ilocal=1,npixbatch
5539 bucketoffsets(ilocal+1)=bucketoffsets(ilocal)+bucketcounts(ilocal)
5541 bucketfill(1:npixbatch)=bucketoffsets(1:npixbatch)
5544 ipix=nint(segments_recv(1,is))
5546 ilocal=ipix-ipixstart+1
5547 idx(bucketfill(ilocal))=is
5548 bucketfill(ilocal)=bucketfill(ilocal)+1
5553 do ipix=ipixstart,ipixend
5555 ilocal=ipix-ipixstart+1
5556 nidx=bucketcounts(ilocal)
5558 profile_local(5)=profile_local(5)+dble(nidx)*dble(nidx)
5560 idx(bucketoffsets(ilocal):bucketoffsets(ilocal+1)-1),nidx)
5561 ixglobal=1+mod(ipix-1,numxi1)
5562 iyglobal=1+(ipix-1)/numxi1
5563 do iseg=bucketoffsets(ilocal),bucketoffsets(ilocal+1)-1
5565 euvthin(ixglobal,iyglobal)=euvthin(ixglobal,iyglobal)+segments_recv(3,is)
5567 euv(ixglobal,iyglobal)=euv(ixglobal,iyglobal)+atten*segments_recv(3,is)
5568 tau(ixglobal,iyglobal)=tau(ixglobal,iyglobal)+max(zero,segments_recv(4,is))
5575 deallocate(segments_send,segments_recv)
5576 if (
allocated(segments))
deallocate(segments)
5580 do iigrid=1,igridstail
5581 if (
allocated(cache(iigrid)%source))
deallocate(cache(iigrid)%source)
5582 if (
allocated(cache(iigrid)%opacity))
deallocate(cache(iigrid)%opacity)
5583 if (
allocated(cache(iigrid)%rface))
deallocate(cache(iigrid)%rface)
5584 if (
allocated(cache(iigrid)%thetaface))
deallocate(cache(iigrid)%thetaface)
5585 if (
allocated(cache(iigrid)%phiface))
deallocate(cache(iigrid)%phiface)
5586 if (
allocated(cache(iigrid)%rface2))
deallocate(cache(iigrid)%rface2)
5587 if (
allocated(cache(iigrid)%theta_cos))
deallocate(cache(iigrid)%theta_cos)
5588 if (
allocated(cache(iigrid)%phi_sin))
deallocate(cache(iigrid)%phi_sin)
5589 if (
allocated(cache(iigrid)%phi_cos))
deallocate(cache(iigrid)%phi_cos)
5592 deallocate(sendcounts,recvcounts,senddispls,recvdispls,roundsendcounts,roundrecvcounts,&
5593 roundsenddispls,roundrecvdispls,ownersegcounts,owneroffsets,bucketcounts,&
5594 bucketoffsets,bucketfill)
5595 allocate(image_reduce(numxi1,numxi2))
5596 call mpi_allreduce(euv,image_reduce,numxi1*numxi2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5598 call mpi_allreduce(tau,image_reduce,numxi1*numxi2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5600 call mpi_allreduce(euvthin,image_reduce,numxi1*numxi2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5601 euvthin=image_reduce
5602 deallocate(image_reduce)
5603 call mpi_allreduce(profile_local,profile_global,5,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5604 call mpi_allreduce(phys_max_local,phys_max_global,2,mpi_double_precision,mpi_max,
icomm,
ierrmpi)
5605 call mpi_allreduce(phys_sum_local,phys_sum_global,2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5606 call mpi_allreduce(sphddafallbacklocal,sphddafallbackglobal,1,mpi_integer,mpi_sum,
icomm,
ierrmpi)
5609 write(*,
'(a,5(es12.5,1x))')
' sph_dda thick profile: ',profile_global
5610 write(*,
'(a,i0)')
' sph_dda thick fallback rays: ',sphddafallbackglobal
5611 write(*,
'(a,4(es12.5,1x))')
' sph_dda thick physics maxj maxk sumjds sumkds: ',&
5612 phys_max_global(1),phys_max_global(2),phys_sum_global(1),phys_sum_global(2)
5614 write(*,
'(a,5(es12.5,1x))')
' sph_intersection thick profile: ',profile_global
5615 write(*,
'(a,4(es12.5,1x))')
' sph_intersection thick physics maxj maxk sumjds sumkds: ',&
5616 phys_max_global(1),phys_max_global(2),phys_sum_global(1),phys_sum_global(2)
5626 double precision,
intent(out) :: xImin1,xImax1,xImin2,xImax2
5628 integer,
parameter :: nsample=5
5629 integer :: iigrid,igrid,ir,it,ip
5630 integer :: ixI^L,ixO^L
5631 double precision,
allocatable :: rface(:),thetaface(:),phiface(:)
5632 double precision :: local_min1,local_max1,local_min2,local_max2
5633 double precision :: sph(1:3),xcent(1:2),wr,wt,wp
5635 local_min1=huge(one)
5636 local_max1=-huge(one)
5637 local_min2=huge(one)
5638 local_max2=-huge(one)
5640 do iigrid=1,igridstail
5641 igrid=igrids(iigrid)
5642 ^d&ixomin^d=ixmlo^d\
5643 ^d&ixomax^d=ixmhi^d\
5644 ^d&iximin^d=ixglo^d\
5645 ^d&iximax^d=ixghi^d\
5648 rface,thetaface,phiface)
5650 wr=dble(ir)/dble(nsample-1)
5651 sph(1)=(one-wr)*rface(ixomin1)+wr*rface(ixomax1+1)
5653 wt=dble(it)/dble(nsample-1)
5654 sph(2)=(one-wt)*thetaface(ixomin2)+wt*thetaface(ixomax2+1)
5656 wp=dble(ip)/dble(nsample-1)
5657 if (ir/=0 .and. ir/=nsample-1 .and. it/=0 .and. it/=nsample-1 .and. &
5658 ip/=0 .and. ip/=nsample-1) cycle
5659 sph(3)=(one-wp)*phiface(ixomin3)+wp*phiface(ixomax3+1)
5661 local_min1=min(local_min1,xcent(1))
5662 local_max1=max(local_max1,xcent(1))
5663 local_min2=min(local_min2,xcent(2))
5664 local_max2=max(local_max2,xcent(2))
5668 deallocate(rface,thetaface,phiface)
5671 call mpi_allreduce(local_min1,ximin1,1,mpi_double_precision,mpi_min,icomm,ierrmpi)
5672 call mpi_allreduce(local_max1,ximax1,1,mpi_double_precision,mpi_max,icomm,ierrmpi)
5673 call mpi_allreduce(local_min2,ximin2,1,mpi_double_precision,mpi_min,icomm,ierrmpi)
5674 call mpi_allreduce(local_max2,ximax2,1,mpi_double_precision,mpi_max,icomm,ierrmpi)
5675 if (ximin1>0.5d0*huge(one) .or. ximax1<-0.5d0*huge(one) .or. &
5676 ximin2>0.5d0*huge(one) .or. ximax2<-0.5d0*huge(one))
then
5677 call mpistop(
"sph_intersection could not determine image bounds")
5682 double precision,
intent(out) :: dxI
5684 select case(trim(dat_resolution_mode))
5690 call mpistop(
"unknown dat_resolution_mode")
5695 double precision,
intent(out) :: dxI
5697 double precision :: refine_factor,dr,dtheta,dphi,rmin,sin_theta_min
5699 refine_factor=dble(2**(refine_max_level-1))
5701 dxi=min(abs(xprobmax1-xprobmin1)/(dble(domain_nx1)*refine_factor),&
5702 abs(xprobmax2-xprobmin2)/(dble(domain_nx2)*refine_factor),&
5703 abs(xprobmax3-xprobmin3)/(dble(domain_nx3)*refine_factor))
5704 else if (coordinate==spherical)
then
5705 rmin=max(smalldouble,min(xprobmin1,xprobmax1))
5706 sin_theta_min=max(smalldouble,min(abs(sin(xprobmin2)),&
5707 abs(sin(xprobmax2))))
5708 dr=abs(xprobmax1-xprobmin1)/(dble(domain_nx1)*refine_factor)
5709 dtheta=abs(xprobmax2-xprobmin2)/(dble(domain_nx2)*refine_factor)
5710 dphi=abs(xprobmax3-xprobmin3)/(dble(domain_nx3)*refine_factor)
5711 dxi=min(dr,rmin*dtheta,rmin*sin_theta_min*dphi)
5713 call mpistop(
"nominal dat resolution needs Cartesian or spherical coordinates")
5716 if (dxi<=zero .or. dxi>half*huge(one))
then
5717 call mpistop(
"could not determine nominal dat-resolution image spacing")
5722 double precision,
intent(out) :: dxI
5724 integer :: iigrid,igrid,ixI^L,ixO^L,ix^D
5725 double precision :: local_min,global_min,dr,ds_theta,ds_phi,rval,theta
5728 do iigrid=1,igridstail
5729 igrid=igrids(iigrid)
5730 ^d&ixomin^d=ixmlo^d\
5731 ^d&ixomax^d=ixmhi^d\
5732 ^d&iximin^d=ixglo^d\
5733 ^d&iximax^d=ixghi^d\
5735 do ix1=ixomin1,ixomax1
5736 do ix2=ixomin2,ixomax2
5737 do ix3=ixomin3,ixomax3
5739 local_min=min(local_min,ps(igrid)%dx(ix^d,1),&
5740 ps(igrid)%dx(ix^d,2),ps(igrid)%dx(ix^d,3))
5741 else if (coordinate==spherical)
then
5742 rval=max(smalldouble,ps(igrid)%x(ix^d,1))
5743 theta=ps(igrid)%x(ix^d,2)
5744 dr=ps(igrid)%dx(ix^d,1)
5745 ds_theta=rval*ps(igrid)%dx(ix^d,2)
5746 ds_phi=rval*max(smalldouble,sin(theta))*ps(igrid)%dx(ix^d,3)
5747 local_min=min(local_min,dr,ds_theta,ds_phi)
5749 call mpistop(
"minimum dat resolution needs Cartesian or spherical coordinates")
5756 call mpi_allreduce(local_min,global_min,1,mpi_double_precision,mpi_min,&
5758 if (global_min<=zero .or. global_min>half*huge(one))
then
5759 call mpistop(
"could not determine minimum dat-resolution image spacing")
5770 integer,
intent(in) :: qunit
5772 character(20),
intent(in) :: datatype
5774 integer :: ix^D,numXI1,numXI2,numWI
5775 double precision :: xImin1,xImax1,xImin2,xImax2,xIcent1,xIcent2,dxI
5776 double precision,
allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:)
5777 double precision,
allocatable :: wI(:,:,:),wIs(:,:,:),EM(:,:),Dpl(:,:),Tau(:,:),EMthin(:,:),WLB(:,:,:)
5778 double precision :: vec_temp1(1:3),vec_temp2(1:3)
5779 double precision :: vec_z(1:3),vec_cor(1:3),xI_cor(1:2)
5780 double precision :: res,LOS_psi,r_max,r_loc
5783 character (30) :: ion
5784 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
5785 double precision :: arcsec,RHESSI_rsl,LASCO_rsl,pixel,R_occult,smallflux
5786 integer :: iigrid,igrid,i,j,numSI,iw
5787 logical :: emit,ray_image_global,has_thick_output
5789 if (coordinate==spherical)
then
5797 if (coordinate==spherical)
then
5801 ximin1=-abs(xprobmax1)
5802 ximin2=-abs(xprobmax1)
5803 ximax1=abs(xprobmax1)
5804 ximax2=abs(xprobmax1)
5809 if (ix1==1) vec_cor(1)=xprobmin1
5810 if (ix1==2) vec_cor(1)=xprobmax1
5812 if (ix2==1) vec_cor(2)=xprobmin2
5813 if (ix2==2) vec_cor(2)=xprobmax2
5815 if (ix3==1) vec_cor(3)=xprobmin3
5816 if (ix3==2) vec_cor(3)=xprobmax3
5819 r_loc=r_loc+(vec_cor(2)-
x_origin(2))**2
5820 r_loc=r_loc+(vec_cor(3)-
x_origin(3))**2
5822 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
5825 r_max=max(r_max,r_loc)
5829 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
5835 ximin1=min(ximin1,xi_cor(1))
5836 ximax1=max(ximax1,xi_cor(1))
5837 ximin2=min(ximin2,xi_cor(2))
5838 ximax2=max(ximax2,xi_cor(2))
5851 xicent1=(ximin1+ximax1)/2.d0
5852 xicent2=(ximin2+ximax2)/2.d0
5860 if (datatype==
'image_euv')
then
5864 else if (datatype==
'image_sxr')
then
5866 dxi=rhessi_rsl*arcsec
5868 else if (datatype==
'image_whitelight')
then
5879 call mpistop(
'Whitelight synthesis: instrument is not supported!')
5881 dxi=lasco_rsl*arcsec
5886 numxi1=8*ceiling((ximax1-xicent1)/dxi/8.d0)
5887 ximin1=xicent1-numxi1*dxi
5888 ximax1=xicent1+numxi1*dxi
5890 numxi2=8*ceiling((ximax2-xicent2)/dxi/8.d0)
5891 ximin2=xicent2-numxi2*dxi
5892 ximax2=xicent2+numxi2*dxi
5894 allocate(xi1(numxi1),xi2(numxi2),dxi1(numxi1),dxi2(numxi2))
5896 xi1(ix1)=ximin1+dxi*(ix1-
half)
5900 xi2(ix2)=ximin2+dxi*(ix2-
half)
5905 if (datatype==
'image_euv' .or. datatype==
'image_sxr')
then
5906 has_thick_output=datatype==
'image_euv' .and. trim(
radiation_transfer)==
'thick' .and. &
5909 if (datatype==
'image_euv')
then
5914 allocate(wi(numxi1,numxi2,numwi),wis(numxi1,numxi2,numwi),em(numxi1,numxi2))
5918 ray_image_global=.false.
5919 if (has_thick_output)
then
5920 allocate(tau(numxi1,numxi2),emthin(numxi1,numxi2))
5924 if (
slab .and. datatype==
'image_euv' .and. &
5926 ray_image_global=.true.
5927 allocate(dpl(numxi1,numxi2))
5936 do iigrid=1,igridstail; igrid=igrids(iigrid);
5939 else if (trim(
ray_method_active) ==
'spherical' .and. datatype ==
'image_euv')
then
5941 ray_image_global=.true.
5947 do iigrid=1,igridstail; igrid=igrids(iigrid);
5951 if (ray_image_global)
then
5952 if (has_thick_output)
then
5954 has_thick_output,tau=tau,euvthin=emthin,&
5955 cap_absorption=.true.)
5962 if (em(ix1,ix2)>smallflux) wis(ix1,ix2,1)=em(ix1,ix2)
5966 if (.not. ray_image_global)
then
5967 numsi=numxi1*numxi2*numwi
5968 call mpi_allreduce(wis,wi,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5976 call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
5977 if (
allocated(tau))
deallocate(tau)
5978 if (
allocated(emthin))
deallocate(emthin)
5979 deallocate(wi,wis,em)
5980 else if (datatype==
'image_whitelight')
then
5982 allocate(wi(numxi1,numxi2,numwi),wis(numxi1,numxi2,numwi),wlb(numxi1,numxi2,numwi))
5986 if (coordinate==spherical)
then
5987 do iigrid=1,igridstail; igrid=igrids(iigrid);
5993 if (wlb(ix1,ix2,1)>smallflux)
then
5994 wis(ix1,ix2,1)=wlb(ix1,ix2,1)
5995 wis(ix1,ix2,2)=wlb(ix1,ix2,2)
5999 numsi=numxi1*numxi2*numwi
6000 call mpi_allreduce(wis,wi,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
6007 call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
6008 deallocate(wi,wis,wlb)
6011 deallocate(xi1,xi2,dxi1,dxi2)
6016 integer,
intent(in) :: igrid,numXI1,numXI2
6017 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
6018 double precision,
intent(in) :: dxI
6020 character(20),
intent(in) :: datatype
6021 double precision,
intent(inout) :: EM(numXI1,numXI2)
6023 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
6024 double precision :: xb^L,xd^D
6025 double precision,
allocatable :: flux(:^D&),opacity(:^D&)
6026 double precision :: res
6027 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
6028 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC
6029 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2)
6032 double precision :: logTe
6033 character (30) :: ion
6034 double precision :: lineCent
6035 double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
6036 double precision :: arcsec,pixel,RHESSI_rsl,area_1AU
6037 double precision :: aa,bb
6039 ^d&ixomin^d=ixmlo^d\
6040 ^d&ixomax^d=ixmhi^d\
6041 ^d&iximin^d=ixglo^d\
6042 ^d&iximax^d=ixghi^d\
6043 ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
6044 ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
6047 arcsec=7.25d5/unit_length
6049 arcsec=7.25d7/unit_length
6052 allocate(flux(ixi^s),opacity(ixi^s))
6053 if (datatype==
'image_euv')
then
6054 if (trim(emission_model)==
'pseudo_current')
then
6056 else if (trim(emission_model)==
'radio_ff')
then
6060 call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
6061 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
6063 call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
6064 pixel=spacersl*arcsec
6065 sigma0=sigma_psf*pixel
6066 else if (datatype==
'image_sxr')
then
6068 call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr)
6069 rhessi_rsl=2.3d0/instrument_resolution_factor
6071 pixel=rhessi_rsl*arcsec
6072 sigma0=sigma_psf*pixel
6077 {
do ix^d=ixomin^d,ixomax^d\}
6079 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi1(^d))/(dxi/2.d0)));
6080 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi2(^d))/(dxi/2.d0)));
6081 ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
6082 if (datatype==
'image_euv')
then
6084 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length*1.d2/dxi/dxi
6086 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length/dxi/dxi
6088 else if (datatype==
'image_sxr')
then
6090 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length**3/area_1au
6092 if (fluxsubc>smalldouble)
then
6094 {
do isubc^d=1,nsubc^d\}
6095 ^d&xsubc(^d)=ps(igrid)%x(ix^dd,^d)-half*ps(igrid)%dx(ix^dd,^d)+(isubc^d-half)*dxsubc^d;
6099 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
6100 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
6101 ixpmin1=max(1,ixp1-3)
6102 ixpmax1=min(ixp1+3,numxi1)
6103 ixpmin2=max(1,ixp2-3)
6104 ixpmax2=min(ixp2+3,numxi2)
6105 do ixp1=ixpmin1,ixpmax1
6106 do ixp2=ixpmin2,ixpmax2
6107 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6108 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6109 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6110 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6111 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
6112 em(ixp1,ixp2)=em(ixp1,ixp2)+fluxsubc*factor
6119 deallocate(flux,opacity)
6123 integer,
intent(in) :: igrid,numXI1,numXI2
6124 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
6125 double precision,
intent(in) :: dxI
6127 character(20),
intent(in) :: datatype
6128 double precision,
intent(inout) :: EM(numXI1,numXI2)
6130 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
6131 double precision,
allocatable :: flux(:^D&),Ne(:^D&),opacity(:^D&)
6132 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
6133 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC,RsubC
6134 double precision :: TBsubC,PBsubC
6135 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2),xSubC_car(1:3)
6136 double precision :: R_thick,dotp,dvolume,R_occult,Rc
6137 double precision :: dxl(1:3),x_sph(1:3),dx_sph(1:3)
6138 double precision :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
6139 logical :: sun_back_side,emit
6142 double precision :: logTe
6143 character (30) :: ion
6144 double precision :: lineCent
6145 double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
6146 double precision :: RHESSI_rsl,area_1AU,arcsec,pixel
6148 ^d&ixomin^d=ixmlo^d;
6149 ^d&ixomax^d=ixmhi^d;
6150 ^d&iximin^d=ixglo^d;
6151 ^d&iximax^d=ixghi^d;
6154 arcsec=7.25d5/unit_length
6156 arcsec=7.25d7/unit_length
6159 allocate(flux(ixi^s),opacity(ixi^s))
6160 if (datatype==
'image_euv')
then
6161 if (trim(emission_model)==
'pseudo_current')
then
6163 else if (trim(emission_model)==
'radio_ff')
then
6167 call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
6168 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
6170 call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
6171 pixel=spacersl*arcsec
6172 sigma0=sigma_psf*pixel
6173 else if (datatype==
'image_sxr')
then
6175 call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr)
6176 rhessi_rsl=2.3d0/instrument_resolution_factor
6178 pixel=rhessi_rsl*arcsec
6179 sigma0=sigma_psf*pixel
6184 r_thick=r_opt_thick*const_rsun/unit_length
6185 {
do ix^d=ixomin^d,ixomax^d\}
6186 x_sph(1:3)=ps(igrid)%x(ix^d,1:3)
6187 dx_sph(1:3)=ps(igrid)%dx(ix^d,1:3)
6189 dxl(2)=x_sph(1)*dx_sph(2)
6190 dxl(3)=x_sph(1)*dsin(x_sph(2))*dx_sph(3)
6195 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
6197 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
6199 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
6201 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
6203 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
6205 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
6210 xsubc(1)=x_sph(1)-half*dx_sph(1)+(isubc1-half)*dx_sph(1)/nsubc1
6212 dxsubc1=dx_sph(1)/nsubc1
6215 xsubc(2)=x_sph(2)-half*dx_sph(2)+(isubc2-half)*dx_sph(2)/nsubc2
6216 dxsubc2=xsubc(1)*dx_sph(2)/nsubc2
6217 dxsubc3=xsubc(1)*dsin(xsubc(2))*dx_sph(3)/nsubc3
6218 dvolume=dxsubc1*dxsubc2*dxsubc3
6219 if (datatype==
'image_euv')
then
6221 fluxsubc=flux(ix^d)*dvolume*unit_length*1.d2/dxi/dxi
6223 fluxsubc=flux(ix^d)*dvolume*unit_length/dxi/dxi
6225 else if (datatype==
'image_sxr')
then
6227 fluxsubc=flux(ix^d)*dvolume*unit_length**3/area_1au
6230 if (fluxsubc>smalldouble)
then
6233 xsubc(3)=x_sph(3)-half*dx_sph(3)+(isubc3-half)*dx_sph(3)/nsubc3
6235 rc=dsqrt(xcent(1)**2+xcent(2)**2)
6240 sun_back_side=.true.
6241 if (dotp<0.d0) sun_back_side=.false.
6243 if (sun_back_side)
then
6245 if (rc>r_thick) emit=.true.
6248 if (xsubc(1)<=r_thick) emit=.false.
6254 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
6255 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
6256 ixpmin1=max(1,ixp1-3)
6257 ixpmax1=min(ixp1+3,numxi1)
6258 ixpmin2=max(1,ixp2-3)
6259 ixpmax2=min(ixp2+3,numxi2)
6260 do ixp1=ixpmin1,ixpmax1
6261 do ixp2=ixpmin2,ixpmax2
6262 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6263 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6264 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6265 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6266 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
6267 em(ixp1,ixp2)=em(ixp1,ixp2)+fluxsubc*factor
6277 deallocate(flux,opacity)
6284 integer,
intent(in) :: igrid,numXI1,numXI2,numWI
6285 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
6286 double precision,
intent(in) :: dxI
6288 character(20),
intent(in) :: datatype
6289 double precision,
intent(inout) :: WLB(numXI1,numXI2,numWI)
6291 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
6292 double precision,
allocatable :: flux(:^D&),Ne(:^D&)
6293 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
6294 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC,RsubC
6295 double precision :: sigma_PSF,sigma0,arcsec,pixel,LASCO_rsl
6296 double precision :: A,B,C,D,Rc,Ne0,TBsubC,PBsubC,factor
6297 double precision :: R_thick,dotp,dvolume,R_occult
6298 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2),xSubC_car(1:3)
6299 double precision :: dxl(1:3),x_sph(1:3),dx_sph(1:3)
6300 double precision :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
6303 ^d&ixomin^d=ixmlo^d;
6304 ^d&ixomax^d=ixmhi^d;
6305 ^d&iximin^d=ixglo^d;
6306 ^d&iximax^d=ixghi^d;
6309 arcsec=7.25d5/unit_length
6311 arcsec=7.25d7/unit_length
6315 if (whitelight_instrument==
'LASCO/C1')
then
6316 lasco_rsl=5.6d0/instrument_resolution_factor
6318 else if (whitelight_instrument==
'LASCO/C2')
then
6319 lasco_rsl=11.4d0/instrument_resolution_factor
6321 else if (whitelight_instrument==
'LASCO/C3')
then
6322 lasco_rsl=56.d0/instrument_resolution_factor
6325 if (r_occultor>1.d0) r_occult=r_occultor
6326 r_occult=r_occult*const_rsun/unit_length
6327 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,ne)
6330 double precision :: nH_dummy(ixI^S)
6331 call eos%get_ne_nH(ixi^l, ixo^l, ps(igrid)%w, ne, nh_dummy)
6334 pixel=lasco_rsl*arcsec
6335 sigma0=sigma_psf*pixel
6338 r_thick=r_opt_thick*const_rsun/unit_length
6339 {
do ix^d=ixomin^d,ixomax^d\}
6340 x_sph(1:3)=ps(igrid)%x(ix^d,1:3)
6341 dx_sph(1:3)=ps(igrid)%dx(ix^d,1:3)
6343 dxl(2)=x_sph(1)*dx_sph(2)
6344 dxl(3)=x_sph(1)*dsin(x_sph(2))*dx_sph(3)
6345 ne0=ne(ix^d)*unit_numberdensity
6350 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
6352 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
6354 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
6356 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
6358 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
6360 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
6365 xsubc(1)=x_sph(1)-half*dx_sph(1)+(isubc1-half)*dx_sph(1)/nsubc1
6367 dxsubc1=dx_sph(1)/nsubc1
6371 xsubc(2)=x_sph(2)-half*dx_sph(2)+(isubc2-half)*dx_sph(2)/nsubc2
6372 dxsubc2=xsubc(1)*dx_sph(2)/nsubc2
6373 dxsubc3=xsubc(1)*dsin(xsubc(2))*dx_sph(3)/nsubc3
6374 dvolume=dxsubc1*dxsubc2*dxsubc3
6377 xsubc(3)=x_sph(3)-half*dx_sph(3)+(isubc3-half)*dx_sph(3)/nsubc3
6379 rc=dsqrt(xcent(1)**2+xcent(2)**2)
6382 if (rc>r_occult)
then
6386 tbsubc=tbsubc*dvolume*unit_length/dxi/dxi
6387 pbsubc=pbsubc*dvolume*unit_length/dxi/dxi
6388 if (tbsubc<1.d-20) emit=.false.
6393 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
6394 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
6395 ixpmin1=max(1,ixp1-3)
6396 ixpmax1=min(ixp1+3,numxi1)
6397 ixpmin2=max(1,ixp2-3)
6398 ixpmax2=min(ixp2+3,numxi2)
6399 do ixp1=ixpmin1,ixpmax1
6400 do ixp2=ixpmin2,ixpmax2
6401 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6402 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6403 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6404 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6405 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
6406 wlb(ixp1,ixp2,1)=wlb(ixp1,ixp2,1)+tbsubc*factor
6407 wlb(ixp1,ixp2,2)=wlb(ixp1,ixp2,2)+pbsubc*factor
6423 double precision,
intent(in) :: Rl
6424 double precision,
intent(inout) :: A,B,C,D
6426 double precision :: sinO,cosO,sinO2,cosO2,tmp
6431 coso=abs(dsqrt(coso2))
6432 tmp=log((1.d0+sino)/coso)
6434 b=-(1.d0-3.d0*sino2-(coso2/sino)*(1.d0+3.d0*sino2)*tmp)/8.d0
6435 c=4.d0/3.d0-coso-coso*coso2/3.d0
6436 d=(5.d0+sino2-(coso2/sino)*(5.d0-sino2)*tmp)/8.d0
6442 double precision,
intent(in) :: Rl,Rin,Ne,A,B,C,D
6443 double precision,
intent(inout) :: fluxTB,fluxPB
6445 double precision :: const,u,Bt,Br,PB,TB,sinchi2
6448 const=1.24878d-25/(1.d0-u/3.d0)
6450 bt=const*(c+u*(d-c))
6451 pb=const*sinchi2*((a+u*(b-a)))
6460 double precision,
intent(in) :: x_sph(1:3)
6461 double precision,
intent(inout) :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
6463 unitv_r(1)=dsin(x_sph(2))*dcos(x_sph(3))
6464 unitv_r(2)=dsin(x_sph(2))*dsin(x_sph(3))
6465 unitv_r(3)=dcos(x_sph(2))
6466 unitv_theta(1)=dcos(x_sph(2))*dcos(x_sph(3))
6467 unitv_theta(2)=dcos(x_sph(2))*dsin(x_sph(3))
6468 unitv_theta(3)=-dsin(x_sph(2))
6469 unitv_phi(1)=-dsin(x_sph(3))
6470 unitv_phi(2)=dcos(x_sph(3))
6475 subroutine output_data(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,datatype)
6479 integer,
intent(in) :: qunit,nXO1,nXO2,nWO
6480 double precision,
intent(in) :: dxO1(nxO1),dxO2(nxO2)
6481 double precision,
intent(in) :: xO1(nXO1),xO2(nxO2)
6482 double precision,
intent(inout) :: wO(nXO1,nXO2,nWO)
6483 character(20),
intent(in) :: datatype
6485 integer :: nPiece,nP1,nP2,nC1,nC2,nWC
6486 integer :: piece_nmax1,piece_nmax2,ix1,ix2,j,ipc,ixc1,ixc2
6487 double precision :: uniform_tol
6488 double precision,
allocatable :: xC(:,:,:,:),wC(:,:,:,:),dxC(:,:,:,:)
6495 if (abs(wo(ix1,ix2,j))<smalldouble) wo(ix1,ix2,j)=zero
6502 if (datatype==
'image_euv' .or. datatype==
'image_sxr')
then
6504 piece_nmax1=block_nx2
6505 piece_nmax2=block_nx3
6507 piece_nmax1=block_nx3
6508 piece_nmax2=block_nx1
6510 piece_nmax1=block_nx1
6511 piece_nmax2=block_nx2
6513 else if (datatype==
'spectrum_euv')
then
6516 piece_nmax2=block_nx1
6518 piece_nmax2=block_nx2
6520 piece_nmax2=block_nx3
6527 loopn1:
do j=piece_nmax1,1,-1
6528 if(mod(nxo1,j)==0)
then
6533 loopn2:
do j=piece_nmax2,1,-1
6534 if(mod(nxo2,j)==0)
then
6547 case(
'EIvtuCCmpi',
'ESvtuCCmpi',
'SIvtuCCmpi',
'WIvtuCCmpi')
6549 allocate(xc(npiece,nc1,nc2,2))
6550 allocate(dxc(npiece,nc1,nc2,2))
6551 allocate(wc(npiece,nc1,nc2,nwo))
6555 ix1=mod(ipc-1,np1)*nc1+ixc1
6556 ix2=floor(1.0*(ipc-1)/np1)*nc2+ixc2
6557 xc(ipc,ixc1,ixc2,1)=xo1(ix1)
6558 xc(ipc,ixc1,ixc2,2)=xo2(ix2)
6559 dxc(ipc,ixc1,ixc2,1)=dxo1(ix1)
6560 dxc(ipc,ixc1,ixc2,2)=dxo2(ix2)
6562 wc(ipc,ixc1,ixc2,j)=wo(ix1,ix2,j)
6569 deallocate(xc,dxc,wc)
6570 case(
'EIvtiCCmpi',
'ESvtiCCmpi',
'SIvtiCCmpi',
'WIvtiCCmpi')
6572 (maxval(abs(dxo1(:)-dxo1(1)))>uniform_tol*max(one,abs(dxo1(1))) .or. &
6573 maxval(abs(dxo2(:)-dxo2(1)))>uniform_tol*max(one,abs(dxo2(1)))))
then
6574 call mpistop(
"vti needs uniform dat-resolution image grids")
6576 call write_image_vticc(qunit,xo1,xo2,dxo1,dxo2,wo,nxo1,nxo2,nwo,nc1,nc2)
6579 call mpistop(
"Error in synthesize emission: Unknown convert_type")
6585 subroutine write_image_vticc(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,nC1,nC2)
6589 integer,
intent(in) :: qunit,nXO1,nXO2,nWO,nC1,nC2
6590 double precision,
intent(in) :: xO1(nXO1),xO2(nxO2)
6591 double precision,
intent(in) :: dxO1(nxO1),dxO2(nxO2)
6592 double precision,
intent(in) :: wO(nXO1,nXO2,nWO)
6594 double precision :: origin(1:3), spacing(1:3)
6595 integer :: wholeExtent(1:6)
6597 integer :: ixC1,ixC2
6601 character (70) :: subname,wname,vname,nameL,nameS
6602 character (len=std_len) :: filename
6603 logical :: sph_datres_no_doppler
6606 origin(1)=xo1(1)-0.5d0*dxo1(1)
6607 origin(2)=xo2(1)-0.5d0*dxo2(1)
6618 inquire(qunit,opened=fileopen)
6619 if(.not.fileopen)
then
6624 write(filename,
'(a,i4.4,a)') trim(
filename_euv),filenr,
".vti"
6626 write(filename,
'(a,i4.4,a)') trim(
filename_sxr),filenr,
".vti"
6632 open(qunit,file=filename,status=
'unknown',form=
'formatted')
6636 write(qunit,
'(a)')
'<?xml version="1.0"?>'
6637 write(qunit,
'(a)',advance=
'no')
'<VTKFile type="ImageData"'
6638 write(qunit,
'(a)')
' version="0.1" byte_order="LittleEndian">'
6639 write(qunit,
'(a,3(1pe14.6),a,6(i10),a,3(1pe14.6),a)')
' <ImageData Origin="',&
6640 origin,
'" WholeExtent="',wholeextent,
'" Spacing="',spacing,
'">'
6642 write(qunit,
'(a)')
'<FieldData>'
6643 write(qunit,
'(2a)')
'<DataArray type="Float32" Name="TIME" ',&
6644 'NumberOfTuples="1" format="ascii">'
6646 write(qunit,
'(a)')
'</DataArray>'
6647 write(qunit,
'(a)')
'</FieldData>'
6649 write(qunit,
'(a,6(i10),a)')
'<Piece Extent="',wholeextent,
'">'
6650 write(qunit,
'(a)')
'<CellData>'
6661 if (trim(
emission_model)==
'pseudo_current' .and. iw==1) vname=
'pseudo_current'
6662 if (trim(
emission_model)==
'radio_ff' .and. iw==1) vname=
'radio_brightness_temperature'
6664 if (iw==2 .and.
dat_resolution .and. (.not. sph_datres_no_doppler) .and. &
6670 ((
dat_resolution .and. ((sph_datres_no_doppler .and. iw==2) .or. &
6671 ((.not. sph_datres_no_doppler) .and. iw==3))) .or. &
6685 vname=
'absorption_fraction'
6696 if (iw==1)
write(vname,
'(a)')
'B'
6697 if (iw==2)
write(vname,
'(a)')
'pB'
6705 write(qunit,
'(a,a,a)')&
6706 '<DataArray type="Float64" Name="',trim(vname),
'" format="ascii">'
6707 write(qunit,
'(200(1pe14.6))') ((wo(ixc1,ixc2,iw),ixc1=1,nxo1),ixc2=1,nxo2)
6708 write(qunit,
'(a)')
'</DataArray>'
6710 write(qunit,
'(a)')
'</CellData>'
6711 write(qunit,
'(a)')
'</Piece>'
6713 write(qunit,
'(a)')
'</ImageData>'
6714 write(qunit,
'(a)')
'</VTKFile>'
6724 integer,
intent(in) :: qunit
6725 integer,
intent(in) :: nPiece,nC1,nC2,nWC
6726 double precision,
intent(in) :: xC(nPiece,nC1,nC2,2),dxC(nPiece,nc1,nc2,2)
6727 double precision,
intent(in) :: wC(nPiece,nC1,nC2,nWC)
6728 character(20),
intent(in) :: datatype
6731 double precision :: xP(nPiece,nC1+1,nC2+1,2)
6734 character (70) :: subname,wname,vname,nameL,nameS
6735 character (len=std_len) :: filename
6736 integer :: ixC1,ixC2,ixP,ix1,ix2,j
6737 integer :: nc,np,icel,VTK_type
6738 logical :: sph_datres_no_doppler
6749 if (ix1<np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1,1,1)-0.5d0*dxc(ixp,ix1,1,1)
6750 if (ix1==np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1-1,1,1)+0.5d0*dxc(ixp,ix1-1,1,1)
6751 if (ix2<np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2,2)-0.5d0*dxc(ixp,1,ix2,2)
6752 if (ix2==np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2-1,2)+0.5d0*dxc(ixp,1,ix2-1,2)
6757 inquire(qunit,opened=fileopen)
6758 if(.not.fileopen)
then
6762 if (datatype==
'image_euv')
then
6763 write(filename,
'(a,i4.4,a)') trim(
filename_euv),filenr,
".vtu"
6764 else if (datatype==
'image_sxr')
then
6765 write(filename,
'(a,i4.4,a)') trim(
filename_sxr),filenr,
".vtu"
6766 else if (datatype==
'image_whitelight')
then
6768 else if (datatype==
'spectrum_euv')
then
6771 open(qunit,file=filename,status=
'unknown',form=
'formatted')
6774 write(qunit,
'(a)')
'<?xml version="1.0"?>'
6775 write(qunit,
'(a)',advance=
'no')
'<VTKFile type="UnstructuredGrid"'
6776 write(qunit,
'(a)')
' version="0.1" byte_order="LittleEndian">'
6777 write(qunit,
'(a)')
'<UnstructuredGrid>'
6778 write(qunit,
'(a)')
'<FieldData>'
6779 write(qunit,
'(2a)')
'<DataArray type="Float32" Name="TIME" ',&
6780 'NumberOfTuples="1" format="ascii">'
6782 write(qunit,
'(a)')
'</DataArray>'
6783 write(qunit,
'(a)')
'</FieldData>'
6785 write(qunit,
'(a,i7,a,i7,a)') &
6786 '<Piece NumberOfPoints="',np,
'" NumberOfCells="',nc,
'">'
6787 write(qunit,
'(a)')
'<CellData>'
6789 if (datatype==
'image_euv')
then
6798 if (trim(
emission_model)==
'pseudo_current') vname=
'pseudo_current'
6799 if (trim(
emission_model)==
'radio_ff') vname=
'radio_brightness_temperature'
6802 if (j==2 .and.
dat_resolution .and. (.not. sph_datres_no_doppler) .and. &
6808 ((
dat_resolution .and. ((sph_datres_no_doppler .and. j==2) .or. &
6809 ((.not. sph_datres_no_doppler) .and. j==3))) .or. &
6823 vname=
'absorption_fraction'
6825 else if (datatype==
'image_sxr')
then
6833 else if (datatype==
'image_whitelight')
then
6834 write(vname,
'(a)')
'whitelight'
6835 else if (datatype==
'spectrum_euv')
then
6842 write(qunit,
'(a,a,a)')&
6843 '<DataArray type="Float64" Name="',trim(vname),
'" format="ascii">'
6844 write(qunit,
'(200(1pe14.6))') ((wc(ixp,ixc1,ixc2,j),ixc1=1,nc1),ixc2=1,nc2)
6845 write(qunit,
'(a)')
'</DataArray>'
6847 write(qunit,
'(a)')
'</CellData>'
6848 write(qunit,
'(a)')
'<Points>'
6849 write(qunit,
'(a)')
'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
6854 write(qunit,
'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
6856 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
6858 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
6862 write(qunit,
'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
6864 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
6866 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
6869 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
6873 write(qunit,
'(a)')
'</DataArray>'
6874 write(qunit,
'(a)')
'</Points>'
6876 write(qunit,
'(a)')
'<Cells>'
6877 write(qunit,
'(a)')
'<DataArray type="Int32" Name="connectivity" format="ascii">'
6880 write(qunit,
'(4(i7))') ix1-1+(ix2-1)*np1,ix1+(ix2-1)*np1,&
6881 ix1-1+ix2*np1,ix1+ix2*np1
6884 write(qunit,
'(a)')
'</DataArray>'
6886 write(qunit,
'(a)')
'<DataArray type="Int32" Name="offsets" format="ascii">'
6888 write(qunit,
'(i7)') icel*(2**2)
6890 write(qunit,
'(a)')
'</DataArray>'
6892 write(qunit,
'(a)')
'<DataArray type="Int32" Name="types" format="ascii">'
6896 write(qunit,
'(i2)') vtk_type
6898 write(qunit,
'(a)')
'</DataArray>'
6899 write(qunit,
'(a)')
'</Cells>'
6900 write(qunit,
'(a)')
'</Piece>'
6902 write(qunit,
'(a)')
'</UnstructuredGrid>'
6903 write(qunit,
'(a)')
'</VTKFile>'
6909 double precision,
intent(in) :: vec1(1:3),vec2(1:3)
6910 double precision,
intent(out) :: res
6912 res=vec1(1)*vec2(1)+vec1(2)*vec2(2)+vec1(3)*vec2(3)
6917 double precision,
intent(in) :: vec_in1(1:3),vec_in2(1:3)
6918 double precision,
intent(out) :: vec_out(1:3)
6920 vec_out(1)=vec_in1(2)*vec_in2(3)-vec_in1(3)*vec_in2(2)
6921 vec_out(2)=vec_in1(3)*vec_in2(1)-vec_in1(1)*vec_in2(3)
6922 vec_out(3)=vec_in1(1)*vec_in2(2)-vec_in1(2)*vec_in2(1)
6928 double precision :: LOS_psi
6929 double precision :: vec_car(1:3),vec_z(1:3),vec_temp1(1:3),vec_temp2(1:3)
6930 double precision :: vec_LOS_sph(1:3),vec_xI1_sph(1:3),vec_xI2_sph(1:3)
6934 vec_los(2)=dpi*los_theta/180.d0
6945 if (los_theta==zero)
then
6947 vec_temp1(2)=dpi/2.d0
6948 vec_temp1(3)=dpi*los_phi/180.d0
6962 los_psi=dpi*image_rotate/180.d0
6963 vec_xi1=vec_temp1*cos(los_psi)-vec_temp2*sin(los_psi)
6964 vec_xi2=vec_temp2*cos(los_psi)+vec_temp1*sin(los_psi)
6975 vec_los_sph(2:3)=vec_los_sph(2:3)*180.d0/dpi
6976 vec_xi1_sph(2:3)=vec_xi1_sph(2:3)*180.d0/dpi
6977 vec_xi2_sph(2:3)=vec_xi2_sph(2:3)*180.d0/dpi
6979 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),
']'
6980 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),
']'
6981 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),
']'
6987 double precision,
intent(in) :: vec_sph(1:3)
6988 double precision,
intent(inout) :: vec_car(1:3)
6990 vec_car(1)=vec_sph(1)*dsin(vec_sph(2))*dcos(vec_sph(3))
6991 vec_car(2)=vec_sph(1)*dsin(vec_sph(2))*dsin(vec_sph(3))
6992 vec_car(3)=vec_sph(1)*dcos(vec_sph(2))
6998 double precision,
intent(in) :: vec_car(1:3)
6999 double precision,
intent(inout) :: vec_sph(1:3)
7001 vec_sph(1)=dsqrt(vec_car(1)**2+vec_car(2)**2+vec_car(3)**2)
7002 vec_sph(2)=dacos(vec_car(3)/vec_sph(1))
7003 vec_sph(3)=atan2(vec_car(2),vec_car(1))
7009 double precision :: LOS_psi
7010 double precision :: vec_z(1:3),vec_temp1(1:3),vec_temp2(1:3)
7013 vec_los(1)=-cos(dpi*los_phi/180.d0)*sin(dpi*los_theta/180.d0)
7014 vec_los(2)=-sin(dpi*los_phi/180.d0)*sin(dpi*los_theta/180.d0)
7015 vec_los(3)=-cos(dpi*los_theta/180.d0)
7021 if (los_theta==zero)
then
7022 vec_xi1(1)=cos(dpi*los_phi/180.d0)
7023 vec_xi1(2)=sin(dpi*los_phi/180.d0)
7031 los_psi=dpi*image_rotate/180.d0
7032 vec_xi1=vec_temp1*cos(los_psi)-vec_temp2*sin(los_psi)
7033 vec_xi2=vec_temp2*cos(los_psi)+vec_temp1*sin(los_psi)
7047 double precision,
intent(in) :: x_3D_sph(1:3)
7048 double precision,
intent(inout) :: x_image(1:2)
7049 double precision :: res,res_origin
7050 double precision :: x_3D(1:3)
7061 double precision,
intent(in) :: x_3D(1:3)
7062 double precision,
intent(inout) :: x_image(1:2)
7063 double precision :: res,res_origin
7067 x_image(1)=res-res_origin
7070 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
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...
type(state), pointer block
Block pointer for using one block and its previous state.
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)