19 double precision ::
t_aia(1:101)
22 double precision ::
f_335(1:101)
35 data t_aia / 4. , 4.05, 4.1, 4.15, 4.2, 4.25, 4.3, 4.35, &
36 4.4, 4.45, 4.5, 4.55, 4.6, 4.65, 4.7, 4.75, &
37 4.8, 4.85, 4.9, 4.95, 5. , 5.05, 5.1, 5.15, &
38 5.2, 5.25, 5.3, 5.35, 5.4, 5.45, 5.5, 5.55, &
39 5.6, 5.65, 5.7, 5.75, 5.8, 5.85, 5.9, 5.95, &
40 6. , 6.05, 6.1, 6.15, 6.2, 6.25, 6.3, 6.35, &
41 6.4, 6.45, 6.5, 6.55, 6.6, 6.65, 6.7, 6.75, &
42 6.8, 6.85, 6.9, 6.95, 7. , 7.05, 7.1, 7.15, &
43 7.2, 7.25, 7.3, 7.35, 7.4, 7.45, 7.5, 7.55, &
44 7.6, 7.65, 7.7, 7.75, 7.8, 7.85, 7.9, 7.95, &
45 8. , 8.05, 8.1, 8.15, 8.2, 8.25, 8.3, 8.35, &
46 8.4, 8.45, 8.5, 8.55, 8.6, 8.65, 8.7, 8.75, &
47 8.8, 8.85, 8.9, 8.95, 9. /
49 data f_94 / 4.25022959
d-37, 4.35880298
d-36, 3.57054296
d-35, 2.18175426
d-34, &
50 8.97592571
d-34, 2.68512961
d-33, 7.49559346
d-33, 2.11603751
d-32, &
51 5.39752853
d-32, 1.02935904
d-31, 1.33822307
d-31, 1.40884290
d-31, &
52 1.54933156
d-31, 2.07543102
d-31, 3.42026227
d-31, 6.31171444
d-31, &
53 1.16559416
d-30, 1.95360497
d-30, 2.77818735
d-30, 3.43552578
d-30, &
54 4.04061803
d-30, 4.75470982
d-30, 5.65553769
d-30, 6.70595782
d-30, &
55 7.80680354
d-30, 8.93247715
d-30, 1.02618156
d-29, 1.25979030
d-29, &
56 1.88526483
d-29, 3.62448572
d-29, 7.50553279
d-29, 1.42337571
d-28, &
57 2.37912813
d-28, 3.55232305
d-28, 4.84985757
d-28, 6.20662827
d-28, &
58 7.66193687
d-28, 9.30403645
d-28, 1.10519802
d-27, 1.25786927
d-27, &
59 1.34362634
d-27, 1.33185242
d-27, 1.22302081
d-27, 1.05677973
d-27, &
60 9.23064720
d-28, 8.78570994
d-28, 8.02397416
d-28, 5.87681142
d-28, &
61 3.82272695
d-28, 3.11492649
d-28, 3.85736090
d-28, 5.98893519
d-28, &
62 9.57553548
d-28, 1.46650267
d-27, 2.10365847
d-27, 2.79406671
d-27, &
63 3.39420087
d-27, 3.71077520
d-27, 3.57296767
d-27, 2.95114380
d-27, &
64 2.02913103
d-27, 1.13361825
d-27, 5.13405629
d-28, 2.01305089
d-28, &
65 8.15781482
d-29, 4.28366817
d-29, 3.08701543
d-29, 2.68693906
d-29, &
66 2.51764203
d-29, 2.41773103
d-29, 2.33996083
d-29, 2.26997246
d-29, &
67 2.20316143
d-29, 2.13810001
d-29, 2.07424438
d-29, 2.01149189
d-29, &
68 1.94980213
d-29, 1.88917920
d-29, 1.82963583
d-29, 1.77116920
d-29, &
69 1.71374392
d-29, 1.65740593
d-29, 1.60214447
d-29, 1.54803205
d-29, &
70 1.49510777
d-29, 1.44346818
d-29, 1.39322305
d-29, 1.34441897
d-29, &
71 1.29713709
d-29, 1.25132618
d-29, 1.20686068
d-29, 1.14226584
d-29, &
72 1.09866413
d-29, 1.05635524
d-29, 1.01532444
d-29, 9.75577134
d-30, &
73 9.37102736
d-30, 8.99873335
d-30, 8.63860172
d-30, 8.29051944
d-30, &
76 data f_131 / 3.18403601
d-37, 3.22254703
d-36, 2.61657920
d-35, &
77 1.59575286
d-34, 6.65779556
d-34, 2.07015132
d-33, &
78 6.05768615
d-33, 1.76074833
d-32, 4.52633001
d-32, &
79 8.57121883
d-32, 1.09184271
d-31, 1.10207963
d-31, &
80 1.11371658
d-31, 1.29105226
d-31, 1.80385897
d-31, &
81 3.27295431
d-31, 8.92002136
d-31, 3.15214579
d-30, &
82 9.73440787
d-30, 2.22709702
d-29, 4.01788984
d-29, &
83 6.27471832
d-29, 8.91764995
d-29, 1.18725647
d-28, &
84 1.52888040
d-28, 2.05082946
d-28, 3.47651873
d-28, &
85 8.80482184
d-28, 2.66533063
d-27, 7.05805149
d-27, &
86 1.46072515
d-26, 2.45282476
d-26, 3.55303726
d-26, &
87 4.59075911
d-26, 5.36503515
d-26, 5.68444094
d-26, &
88 5.47222296
d-26, 4.81119761
d-26, 3.85959059
d-26, &
89 2.80383406
d-26, 1.83977650
d-26, 1.11182849
d-26, &
90 6.50748885
d-27, 3.96843481
d-27, 2.61876319
d-27, &
91 1.85525324
d-27, 1.39717024
d-27, 1.11504283
d-27, &
92 9.38169611
d-28, 8.24801234
d-28, 7.43331919
d-28, &
93 6.74537063
d-28, 6.14495760
d-28, 5.70805277
d-28, &
94 5.61219786
d-28, 6.31981777
d-28, 9.19747307
d-28, &
95 1.76795732
d-27, 3.77985446
d-27, 7.43166191
d-27, &
96 1.19785603
d-26, 1.48234676
d-26, 1.36673114
d-26, &
97 9.61047146
d-27, 5.61209353
d-27, 3.04779780
d-27, &
98 1.69378976
d-27, 1.02113491
d-27, 6.82223774
d-28, &
99 5.02099099
d-28, 3.99377760
d-28, 3.36279037
d-28, &
100 2.94767378
d-28, 2.65740865
d-28, 2.44396277
d-28, &
101 2.28003967
d-28, 2.14941419
d-28, 2.04178995
d-28, &
102 1.95031045
d-28, 1.87011994
d-28, 1.79777869
d-28, &
103 1.73093957
d-28, 1.66795789
d-28, 1.60785455
d-28, &
104 1.55002399
d-28, 1.49418229
d-28, 1.44022426
d-28, &
105 1.38807103
d-28, 1.33772767
d-28, 1.28908404
d-28, &
106 1.24196208
d-28, 1.17437501
d-28, 1.12854330
d-28, &
107 1.08410498
d-28, 1.04112003
d-28, 9.99529904
d-29, &
108 9.59358806
d-29, 9.20512291
d-29, 8.83009123
d-29, &
109 8.46817043
d-29, 8.11921928
d-29 /
111 data f_171 / 2.98015581
d-42, 1.24696230
d-40, 3.37614652
d-39, 5.64103034
d-38, &
112 5.20550266
d-37, 2.77785939
d-36, 1.16283616
d-35, 6.50007689
d-35, &
113 9.96177399
d-34, 1.89586076
d-32, 2.10982799
d-31, 1.36946479
d-30, &
114 6.27396553
d-30, 2.29955134
d-29, 7.13430211
d-29, 1.91024282
d-28, &
115 4.35358848
d-28, 7.94807808
d-28, 1.07431875
d-27, 1.08399488
d-27, &
116 9.16212938
d-28, 7.34715770
d-28, 6.59246382
d-28, 9.13541375
d-28, &
117 2.05939035
d-27, 5.08206555
d-27, 1.10148083
d-26, 2.01884662
d-26, &
118 3.13578384
d-26, 4.14367719
d-26, 5.36067711
d-26, 8.74170213
d-26, &
119 1.64161233
d-25, 2.94587860
d-25, 4.76298332
d-25, 6.91765639
d-25, &
120 9.08825111
d-25, 1.08496183
d-24, 1.17440114
d-24, 1.13943939
d-24, &
121 9.71696981
d-25, 7.09593688
d-25, 4.31376399
d-25, 2.12708486
d-25, &
122 8.47429567
d-26, 3.17608104
d-26, 1.95898842
d-26, 1.98064242
d-26, &
123 1.67706555
d-26, 8.99126003
d-27, 3.29773878
d-27, 1.28896127
d-27, &
124 8.51169698
d-28, 7.53520167
d-28, 6.18268143
d-28, 4.30034650
d-28, &
125 2.78152409
d-28, 1.95437088
d-28, 1.65896278
d-28, 1.68740181
d-28, &
126 1.76054383
d-28, 1.63978419
d-28, 1.32880591
d-28, 1.00833205
d-28, &
127 7.82252806
d-29, 6.36181741
d-29, 5.34633869
d-29, 4.58013864
d-29, &
128 3.97833422
d-29, 3.49414760
d-29, 3.09790940
d-29, 2.76786227
d-29, &
129 2.48806269
d-29, 2.24823367
d-29, 2.04016653
d-29, 1.85977413
d-29, &
130 1.70367499
d-29, 1.56966125
d-29, 1.45570643
d-29, 1.35964565
d-29, &
131 1.27879263
d-29, 1.21016980
d-29, 1.15132499
d-29, 1.09959628
d-29, &
132 1.05307482
d-29, 1.01040261
d-29, 9.70657096
d-30, 9.33214234
d-30, &
133 8.97689427
d-30, 8.63761192
d-30, 8.31149879
d-30, 7.85162401
d-30, &
134 7.53828281
d-30, 7.23559452
d-30, 6.94341530
d-30, 6.66137038
d-30, &
135 6.38929156
d-30, 6.12669083
d-30, 5.87346434
d-30, 5.62943622
d-30, &
138 data f_193 / 6.40066486
d-32, 4.92737300
d-31, 2.95342934
d-30, 1.28061594
d-29, &
139 3.47747667
d-29, 5.88554792
d-29, 7.72171179
d-29, 9.75609282
d-29, &
140 1.34318963
d-28, 1.96252638
d-28, 2.70163878
d-28, 3.63192965
d-28, &
141 5.28087341
d-28, 8.37821446
d-28, 1.39089159
d-27, 2.31749718
d-27, &
142 3.77510689
d-27, 5.85198594
d-27, 8.26021568
d-27, 1.04870405
d-26, &
143 1.25209374
d-26, 1.47406787
d-26, 1.77174067
d-26, 2.24098537
d-26, &
144 3.05926105
d-26, 4.50018853
d-26, 6.84720216
d-26, 1.00595861
d-25, &
145 1.30759222
d-25, 1.36481773
d-25, 1.15943558
d-25, 1.01467304
d-25, &
146 1.04092532
d-25, 1.15071251
d-25, 1.27416033
d-25, 1.38463476
d-25, &
147 1.47882726
d-25, 1.57041238
d-25, 1.69786224
d-25, 1.94970397
d-25, &
148 2.50332918
d-25, 3.58321431
d-25, 5.18061550
d-25, 6.60405549
d-25, &
149 6.64085365
d-25, 4.83825816
d-25, 2.40545020
d-25, 8.59534098
d-26, &
150 2.90920638
d-26, 1.33204845
d-26, 9.03933926
d-27, 7.78910836
d-27, &
151 7.29342321
d-27, 7.40267022
d-27, 8.05279981
d-27, 8.13829291
d-27, &
152 6.92634262
d-27, 5.12521880
d-27, 3.59527615
d-27, 2.69617560
d-27, &
153 2.84432713
d-27, 5.06697306
d-27, 1.01281903
d-26, 1.63526978
d-26, &
154 2.06759342
d-26, 2.19482312
d-26, 2.10050611
d-26, 1.89837248
d-26, &
155 1.66347131
d-26, 1.43071097
d-26, 1.21518419
d-26, 1.02078343
d-26, &
156 8.46936184
d-27, 6.93015742
d-27, 5.56973237
d-27, 4.38951754
d-27, &
157 3.38456457
d-27, 2.55309556
d-27, 1.88904224
d-27, 1.38057546
d-27, &
158 1.00718330
d-27, 7.43581116
d-28, 5.63562931
d-28, 4.43359435
d-28, &
159 3.63923535
d-28, 3.11248143
d-28, 2.75586846
d-28, 2.50672237
d-28, &
160 2.32419348
d-28, 2.18325682
d-28, 2.06834486
d-28, 1.93497044
d-28, &
161 1.84540751
d-28, 1.76356504
d-28, 1.68741425
d-28, 1.61566157
d-28, &
162 1.54754523
d-28, 1.48249410
d-28, 1.42020176
d-28, 1.36045230
d-28, &
165 data f_211 / 4.74439912
d-42, 1.95251522
d-40, 5.19700194
d-39, 8.53120166
d-38, &
166 7.72745727
d-37, 4.04158559
d-36, 1.64853511
d-35, 8.56295439
d-35, &
167 1.17529722
d-33, 2.16867729
d-32, 2.40472264
d-31, 1.56418133
d-30, &
168 7.20032889
d-30, 2.65838271
d-29, 8.33196904
d-29, 2.26128236
d-28, &
169 5.24295811
d-28, 9.77791121
d-28, 1.35913489
d-27, 1.43957785
d-27, &
170 1.37591544
d-27, 1.49029886
d-27, 2.06183401
d-27, 3.31440622
d-27, &
171 5.42497318
d-27, 8.41100374
d-27, 1.17941366
d-26, 1.49269794
d-26, &
172 1.71506074
d-26, 1.71266353
d-26, 1.51434781
d-26, 1.36766622
d-26, &
173 1.33483562
d-26, 1.36834518
d-26, 1.45829002
d-26, 1.62575306
d-26, &
174 1.88773347
d-26, 2.22026986
d-26, 2.54930499
d-26, 2.80758138
d-26, &
175 3.06176409
d-26, 3.62799792
d-26, 5.13226109
d-26, 8.46260744
d-26, &
176 1.38486586
d-25, 1.86192535
d-25, 1.78007934
d-25, 1.16548409
d-25, &
177 5.89293257
d-26, 2.69952884
d-26, 1.24891081
d-26, 6.41273176
d-27, &
178 4.08282914
d-27, 3.26463328
d-27, 2.76230280
d-27, 2.08986882
d-27, &
179 1.37658470
d-27, 8.48489381
d-28, 5.19304217
d-28, 3.19312514
d-28, &
180 2.02968197
d-28, 1.50171666
d-28, 1.39164218
d-28, 1.42448821
d-28, &
181 1.41714519
d-28, 1.33341059
d-28, 1.20759270
d-28, 1.07259692
d-28, &
182 9.44895400
d-29, 8.29030041
d-29, 7.25440631
d-29, 6.33479483
d-29, &
183 5.51563757
d-29, 4.79002469
d-29, 4.14990482
d-29, 3.59384972
d-29, &
184 3.12010860
d-29, 2.72624742
d-29, 2.40734791
d-29, 2.15543565
d-29, &
185 1.95921688
d-29, 1.80682882
d-29, 1.68695662
d-29, 1.59020936
d-29, &
186 1.50940886
d-29, 1.43956179
d-29, 1.37731622
d-29, 1.32049043
d-29, &
187 1.26771875
d-29, 1.21803879
d-29, 1.17074716
d-29, 1.10507836
d-29, &
188 1.06022834
d-29, 1.01703080
d-29, 9.75436986
d-30, 9.35349257
d-30, &
189 8.96744546
d-30, 8.59527489
d-30, 8.23678940
d-30, 7.89144480
d-30, &
192 data f_304 / 3.62695850
d-32, 2.79969087
d-31, 1.68340584
d-30, 7.32681440
d-30, &
193 1.99967770
d-29, 3.41296785
d-29, 4.55409104
d-29, 5.94994635
d-29, &
194 8.59864963
d-29, 1.39787633
d-28, 3.17701965
d-28, 1.14474920
d-27, &
195 4.44845958
d-27, 1.54785841
d-26, 4.70265345
d-26, 1.24524365
d-25, &
196 2.81535352
d-25, 5.10093666
d-25, 6.83545307
d-25, 6.82110329
d-25, &
197 5.66886188
d-25, 4.36205513
d-25, 3.29265688
d-25, 2.49802368
d-25, &
198 1.92527113
d-25, 1.51058572
d-25, 1.20596047
d-25, 9.76884267
d-26, &
199 7.89979266
d-26, 6.18224289
d-26, 4.67298332
d-26, 3.57934505
d-26, &
200 2.84535785
d-26, 2.32853022
d-26, 1.95228514
d-26, 1.67880071
d-26, &
201 1.47608785
d-26, 1.32199691
d-26, 1.20070960
d-26, 1.09378177
d-26, &
202 1.00031730
d-26, 9.62434001
d-27, 1.05063954
d-26, 1.27267143
d-26, &
203 1.45923057
d-26, 1.36746707
d-26, 1.03466970
d-26, 6.97647829
d-27, &
204 4.63141039
d-27, 3.19031994
d-27, 2.33373613
d-27, 1.81589079
d-27, &
205 1.48446917
d-27, 1.26611478
d-27, 1.12617468
d-27, 1.03625148
d-27, &
206 9.61400595
d-28, 8.79016231
d-28, 7.82612130
d-28, 6.73762960
d-28, &
207 5.59717956
d-28, 4.53010243
d-28, 3.65712196
d-28, 3.00958686
d-28, &
208 2.54011502
d-28, 2.18102277
d-28, 1.88736437
d-28, 1.63817539
d-28, &
209 1.42283147
d-28, 1.23631916
d-28, 1.07526003
d-28, 9.36797928
d-29, &
210 8.18565660
d-29, 7.18152734
d-29, 6.32523238
d-29, 5.59513985
d-29, &
211 4.96614048
d-29, 4.42518826
d-29, 3.95487628
d-29, 3.54690294
d-29, &
212 3.18953930
d-29, 2.87720933
d-29, 2.60186750
d-29, 2.36011522
d-29, &
213 2.14717806
d-29, 1.95905217
d-29, 1.79287981
d-29, 1.64562262
d-29, &
214 1.51489425
d-29, 1.39876064
d-29, 1.29496850
d-29, 1.18665438
d-29, &
215 1.10240474
d-29, 1.02643099
d-29, 9.57780996
d-30, 8.95465151
d-30, &
216 8.38950190
d-30, 7.87283711
d-30, 7.40136507
d-30, 6.96804279
d-30, &
219 data f_335 / 2.46882661
d-32, 1.89476632
d-31, 1.13216502
d-30, 4.89532008
d-30, &
220 1.32745970
d-29, 2.25390335
d-29, 3.00511672
d-29, 3.96035934
d-29, &
221 5.77977656
d-29, 8.58600736
d-29, 1.14083000
d-28, 1.48644411
d-28, &
222 2.15788823
d-28, 3.51628877
d-28, 6.12200698
d-28, 1.08184987
d-27, &
223 1.85590697
d-27, 2.91679107
d-27, 3.94405396
d-27, 4.63610680
d-27, &
224 5.13824456
d-27, 5.66602209
d-27, 6.30009232
d-27, 7.03422868
d-27, &
225 7.77973918
d-27, 8.32371831
d-27, 8.56724316
d-27, 8.62601374
d-27, &
226 8.13308844
d-27, 6.53188216
d-27, 4.55197029
d-27, 3.57590087
d-27, &
227 3.59571707
d-27, 4.03502770
d-27, 4.54366411
d-27, 4.96914990
d-27, &
228 5.24601170
d-27, 5.39979250
d-27, 5.43023669
d-27, 5.26235042
d-27, &
229 4.91585495
d-27, 4.52628362
d-27, 4.13385020
d-27, 3.67538967
d-27, &
230 3.39939742
d-27, 3.81284533
d-27, 5.02332701
d-27, 6.19438602
d-27, &
231 6.49613071
d-27, 6.04010475
d-27, 5.24664275
d-27, 4.37225997
d-27, &
232 3.52957182
d-27, 2.76212276
d-27, 2.08473158
d-27, 1.50850518
d-27, &
233 1.04602472
d-27, 7.13091243
d-28, 5.34289645
d-28, 5.21079581
d-28, &
234 6.22246365
d-28, 6.99555864
d-28, 6.29665489
d-28, 4.45077026
d-28, &
235 2.67046793
d-28, 1.52774686
d-28, 9.18061770
d-29, 6.09116074
d-29, &
236 4.48562572
d-29, 3.59463696
d-29, 3.05820218
d-29, 2.70766652
d-29, &
237 2.46144034
d-29, 2.27758450
d-29, 2.13331183
d-29, 2.01537836
d-29, &
238 1.91566180
d-29, 1.82893912
d-29, 1.75167748
d-29, 1.68136168
d-29, &
239 1.61615595
d-29, 1.55481846
d-29, 1.49643236
d-29, 1.44046656
d-29, &
240 1.38657085
d-29, 1.33459068
d-29, 1.28447380
d-29, 1.23615682
d-29, &
241 1.18963296
d-29, 1.14478976
d-29, 1.10146637
d-29, 1.04039479
d-29, &
242 9.98611410
d-30, 9.58205147
d-30, 9.19202009
d-30, 8.81551313
d-30, &
243 8.45252127
d-30, 8.10224764
d-30, 7.76469090
d-30, 7.43954323
d-30, &
249 data t_iris / 4. , 4.1 , 4.2 , 4.3 , 4.40000001, &
250 4.50000001, 4.60000001, 4.70000001, 4.80000001, 4.90000001, &
251 5.00000001, 5.10000002, 5.20000002, 5.30000002, 5.40000002, &
252 5.50000002, 5.60000002, 5.70000003, 5.80000003, 5.90000003, &
253 6.00000003, 6.10000003, 6.20000003, 6.30000003, 6.40000004, &
254 6.50000004, 6.60000004, 6.70000004, 6.80000004, 6.90000004, &
255 7.00000004, 7.10000005, 7.20000005, 7.30000005, 7.40000005, &
256 7.50000005, 7.60000005, 7.70000006, 7.80000006, 7.90000006, &
259 data f_1354 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
260 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
261 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
262 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
263 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
264 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 1.09503647
d-39, &
265 5.47214550
d-36, 2.42433983
d-33, 2.75295034
d-31, 1.21929718
d-29, &
266 2.48392125
d-28, 2.33268145
d-27, 8.68623633
d-27, 1.00166284
d-26, &
267 3.63126633
d-27, 7.45174807
d-28, 1.38224064
d-28, 2.69270994
d-29, &
268 5.53314977
d-30, 1.15313092
d-30, 2.34195788
d-31, 4.48242942
d-32, &
274 data t_eis1 / 1.99526231
d+05, 2.23872114
d+05, 2.51188643
d+05, 2.81838293
d+05, &
275 3.16227766
d+05, 3.54813389
d+05, 3.98107171
d+05, 4.46683592
d+05, &
276 5.01187234
d+05, 5.62341325
d+05, 6.30957344
d+05, 7.07945784
d+05, &
277 7.94328235
d+05, 8.91250938
d+05, 1.00000000
d+06, 1.12201845
d+06, &
278 1.25892541
d+06, 1.41253754
d+06, 1.58489319
d+06, 1.77827941
d+06, &
279 1.99526231
d+06, 2.23872114
d+06, 2.51188643
d+06, 2.81838293
d+06, &
280 3.16227766
d+06, 3.54813389
d+06, 3.98107171
d+06, 4.46683592
d+06, &
281 5.01187234
d+06, 5.62341325
d+06, 6.30957344
d+06, 7.07945784
d+06, &
282 7.94328235
d+06, 8.91250938
d+06, 1.00000000
d+07, 1.12201845
d+07, &
283 1.25892541
d+07, 1.41253754
d+07, 1.58489319
d+07, 1.77827941
d+07, &
284 1.99526231
d+07, 2.23872114
d+07, 2.51188643
d+07, 2.81838293
d+07, &
285 3.16227766
d+07, 3.54813389
d+07, 3.98107171
d+07, 4.46683592
d+07, &
286 5.01187234
d+07, 5.62341325
d+07, 6.30957344
d+07, 7.07945784
d+07, &
287 7.94328235
d+07, 8.91250938
d+07, 1.00000000
d+08, 1.12201845
d+08, &
288 1.25892541
d+08, 1.41253754
d+08, 1.58489319
d+08, 1.77827941
d+08 /
290 data t_eis2 / 1.99526231
d+06, 2.23872114
d+06, 2.51188643
d+06, 2.81838293
d+06, &
291 3.16227766
d+06, 3.54813389
d+06, 3.98107171
d+06, 4.46683592
d+06, &
292 5.01187234
d+06, 5.62341325
d+06, 6.30957344
d+06, 7.07945784
d+06, &
293 7.94328235
d+06, 8.91250938
d+06, 1.00000000
d+07, 1.12201845
d+07, &
294 1.25892541
d+07, 1.41253754
d+07, 1.58489319
d+07, 1.77827941
d+07, &
295 1.99526231
d+07, 2.23872114
d+07, 2.51188643
d+07, 2.81838293
d+07, &
296 3.16227766
d+07, 3.54813389
d+07, 3.98107171
d+07, 4.46683592
d+07, &
297 5.01187234
d+07, 5.62341325
d+07, 6.30957344
d+07, 7.07945784
d+07, &
298 7.94328235
d+07, 8.91250938
d+07, 1.00000000
d+08, 1.12201845
d+08, &
299 1.25892541
d+08, 1.41253754
d+08, 1.58489319
d+08, 1.77827941
d+08, &
300 1.99526231
d+08, 2.23872114
d+08, 2.51188643
d+08, 2.81838293
d+08, &
301 3.16227766
d+08, 3.54813389
d+08, 3.98107171
d+08, 4.46683592
d+08, &
302 5.01187234
d+08, 5.62341325
d+08, 6.30957344
d+08, 7.07945784
d+08, &
303 7.94328235
d+08, 8.91250938
d+08, 1.00000000
d+09, 1.12201845
d+09, &
304 1.25892541
d+09, 1.41253754
d+09, 1.58489319
d+09, 1.77827941
d+09 /
306 data f_263 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
307 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
308 0.00000000
d+00, 4.46454917
d-45, 3.26774829
d-42, 1.25292566
d-39, &
309 2.66922338
d-37, 3.28497742
d-35, 2.38677554
d-33, 1.03937729
d-31, &
310 2.75075687
d-30, 4.47961733
d-29, 4.46729177
d-28, 2.64862689
d-27, &
311 8.90863800
d-27, 1.72437548
d-26, 2.22217752
d-26, 2.27999477
d-26, &
312 2.08264363
d-26, 1.78226687
d-26, 1.45821699
d-26, 1.14675379
d-26, &
313 8.63082492
d-27, 6.15925429
d-27, 4.11252514
d-27, 2.51530564
d-27, &
314 1.37090986
d-27, 6.42443134
d-28, 2.48392636
d-28, 7.59187874
d-29, &
315 1.77852938
d-29, 3.23945221
d-30, 4.90533903
d-31, 6.75458158
d-32, &
316 9.06878868
d-33, 1.23927474
d-33, 1.75769395
d-34, 2.60710914
d-35, &
317 4.04318030
d-36, 6.53500581
d-37, 1.09365022
d-37, 1.88383322
d-38, &
318 3.31425233
d-39, 5.90964084
d-40, 1.06147549
d-40, 1.90706170
d-41, &
319 3.41331584
d-42, 6.07310718
d-43, 1.07364738
d-43, 1.89085498
d-44, &
320 3.32598922
d-45, 5.87125640
d-46, 0.00000000
d+00, 0.00000000
d+00 /
322 data f_264 / 0.00000000
d+00, 2.81670057
d-46, 1.28007268
d-43, 2.54586603
d-41, &
323 2.67887256
d-39, 1.68413285
d-37, 6.85702304
d-36, 1.91797284
d-34, &
324 3.84675839
d-33, 5.69939170
d-32, 6.36224608
d-31, 5.39176489
d-30, &
325 3.45478458
d-29, 1.64848693
d-28, 5.71476364
d-28, 1.39909997
d-27, &
326 2.37743056
d-27, 2.86712530
d-27, 2.65206348
d-27, 2.07175767
d-27, &
327 1.47866767
d-27, 1.01087374
d-27, 6.79605811
d-28, 4.54746770
d-28, &
328 3.04351751
d-28, 2.03639149
d-28, 1.35940991
d-28, 9.01451939
d-29, &
329 5.91289972
d-29, 3.81821178
d-29, 2.41434696
d-29, 1.48871220
d-29, &
330 8.93362094
d-30, 5.21097445
d-30, 2.95964719
d-30, 1.64278748
d-30, &
331 8.95571660
d-31, 4.82096011
d-31, 2.57390991
d-31, 1.36821781
d-31, &
332 7.27136350
d-32, 3.87019426
d-32, 2.06883430
d-32, 1.11228884
d-32, &
333 6.01883313
d-33, 3.27790676
d-33, 1.79805012
d-33, 9.93085346
d-34, &
334 5.52139556
d-34, 3.08881387
d-34, 1.73890315
d-34, 9.84434964
d-35, &
335 5.60603378
d-35, 3.20626492
d-35, 1.84111068
d-35, 0.00000000
d+00, &
336 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00 /
338 data f_192 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 4.35772105
d-44, &
339 1.26162319
d-41, 1.97471205
d-39, 1.83409019
d-37, 1.08206288
d-35, &
340 4.27914363
d-34, 1.17943846
d-32, 2.32565755
d-31, 3.33087991
d-30, &
341 3.47013260
d-29, 2.60375866
d-28, 1.37737127
d-27, 5.01053913
d-27, &
342 1.23479810
d-26, 2.11310542
d-26, 2.71831513
d-26, 2.89851163
d-26, &
343 2.77312376
d-26, 2.50025229
d-26, 2.18323661
d-26, 1.86980322
d-26, &
344 1.58035034
d-26, 1.31985651
d-26, 1.08733133
d-26, 8.81804906
d-27, &
345 7.00417973
d-27, 5.43356567
d-27, 4.09857884
d-27, 2.99651764
d-27, &
346 2.11902962
d-27, 1.45014127
d-27, 9.62291023
d-28, 6.21548647
d-28, &
347 3.92807578
d-28, 2.44230375
d-28, 1.50167782
d-28, 9.17611405
d-29, &
348 5.58707641
d-29, 3.40570915
d-29, 2.08030862
d-29, 1.27588676
d-29, &
349 7.86535588
d-30, 4.87646151
d-30, 3.03888897
d-30, 1.90578649
d-30, &
350 1.20195947
d-30, 7.61955060
d-31, 4.85602199
d-31, 3.11049969
d-31, &
351 2.00087065
d-31, 1.29223740
d-31, 8.37422008
d-32, 0.00000000
d+00, &
352 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00 /
354 data f_255 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 1.76014287
d-44, &
355 5.07057938
d-42, 7.90473970
d-40, 7.31852999
d-38, 4.30709255
d-36, &
356 1.70009061
d-34, 4.67925160
d-33, 9.21703546
d-32, 1.31918676
d-30, &
357 1.37393161
d-29, 1.03102379
d-28, 5.45694018
d-28, 1.98699648
d-27, &
358 4.90346776
d-27, 8.40524725
d-27, 1.08321456
d-26, 1.15714525
d-26, &
359 1.10905152
d-26, 1.00155023
d-26, 8.75799161
d-27, 7.50935839
d-27, &
360 6.35253533
d-27, 5.30919268
d-27, 4.37669455
d-27, 3.55185164
d-27, &
361 2.82347055
d-27, 2.19257595
d-27, 1.65589541
d-27, 1.21224987
d-27, &
362 8.58395132
d-28, 5.88163935
d-28, 3.90721447
d-28, 2.52593407
d-28, &
363 1.59739995
d-28, 9.93802874
d-29, 6.11343388
d-29, 3.73711135
d-29, &
364 2.27618743
d-29, 1.38793199
d-29, 8.48060787
d-30, 5.20305940
d-30, &
365 3.20867365
d-30, 1.99011277
d-30, 1.24064551
d-30, 7.78310544
d-31, &
366 4.91013681
d-31, 3.11338381
d-31, 1.98451675
d-31, 1.27135460
d-31, &
367 8.17917486
d-32, 5.28280497
d-32, 3.42357159
d-32, 0.00000000
d+00, &
368 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00 /
373 integer,
intent(in) :: ixI^L, ixO^L
374 double precision,
intent(in) :: w(ixI^S,nw)
375 double precision,
intent(in) :: x(ixI^S,1:ndim)
376 double precision,
intent(out):: res(ixI^S)
384 integer,
intent(in) :: ixI^L, ixO^L
385 double precision,
intent(in) :: w(ixI^S, nw)
386 double precision,
intent(in) :: x(ixI^S, 1:ndim)
387 double precision,
intent(out):: val1(ixI^S), val2(ixI^S)
393 procedure(
get_subr1),
pointer,
nopass :: get_rho => null()
394 procedure(
get_subr1),
pointer,
nopass :: get_pthermal => null()
395 procedure(
get_subr1),
pointer,
nopass :: get_var_rfactor => null()
406 double precision,
allocatable :: source(:^d&)
407 double precision,
allocatable :: opacity(:^d&)
408 double precision,
allocatable :: sourcev(:^d&)
409 double precision,
allocatable :: xface1(:),xface2(:),xface3(:)
410 double precision,
allocatable :: rface(:),thetaface(:),phiface(:)
411 double precision,
allocatable :: rface2(:),theta_cos(:),phi_sin(:),phi_cos(:)
412 double precision :: box_min(1:3)=0.d0
413 double precision :: box_max(1:3)=0.d0
418 logical :: has_pixels=.false.
430 character(len=*),
intent(in) :: datatype
433 call mpistop(
"bad radiation_transfer")
438 call mpistop(
"dat_resolution_mode must be nominal or minimum")
454 case(
'cart',
'cart_dda')
456 case(
'spherical',
'sph_intersection')
458 case(
'sph_dda',
'spherical_dda')
468 call mpistop(
"bad emission_model")
472 call mpistop(
"tau and absorption-fraction output need thick transfer")
476 call mpistop(
"radsyn_pixel_batch must be positive")
479 call mpistop(
"radsyn_segment_batch_factor must be non-negative")
482 call mpistop(
"radsyn_segment_memory_mb must be positive")
485 call mpistop(
"radsyn_segment_comm_factor must be positive")
490 call mpistop(
"instrument_postprocess currently needs dat-resolution EUV images")
493 call mpistop(
"instrument_postprocess is not yet supported for spherical rays")
496 call mpistop(
"instrument_postprocess currently supports only EUV AIA or radio_ff images")
499 call mpistop(
"radio_ff instrument_postprocess needs radio_beam_fwhm > 0 arcsec")
507 if (datatype /=
'image_euv' .and. datatype /=
'spectrum_euv')
then
508 call mpistop(
"emission_model=euv_aia is only valid for EUV synthesis")
511 if (datatype /=
'image_whitelight')
then
512 call mpistop(
"emission_model=white_light is only valid for white-light synthesis")
515 if (datatype /=
'image_euv')
then
516 call mpistop(
"emission_model=radio_ff currently reuses EUV-image convert types")
519 call mpistop(
"emission_model=radio_ff needs radio_frequency > 0")
521 case(
'pseudo_current')
522 if (datatype /=
'image_euv')
then
523 call mpistop(
"emission_model=pseudo_current is only valid for EUV-image convert types")
526 call mpistop(
"emission_model=pseudo_current currently supports only thin transfer")
531 if (datatype /=
'image_euv' .or. .not.
slab)
then
532 call mpistop(
"ray_method=cart needs Cartesian EUV slab images")
537 call mpistop(
"ray_method=spherical currently needs 3D spherical grids")
540 call mpistop(
"ray_method=spherical currently needs 3D spherical grids")
545 call mpistop(
"bad ray_method=spherical mode")
548 call mpistop(
"ray_method=spherical currently supports only EUV AIA emission")
550 if (xprobmin2<=1.
d-10 .or. xprobmax2>=dpi-1.
d-10)
then
551 call mpistop(
"ray_method=spherical does not support polar-axis crossing domains")
553 if (xprobmax3<=xprobmin3 .or. xprobmax3-xprobmin3>=2.d0*dpi-1.
d-10)
then
554 call mpistop(
"ray_method=spherical does not support phi-wrapping domains")
559 if (datatype /=
'image_euv')
then
560 call mpistop(
"thick transfer is only defined for EUV images")
565 if (.not.
slab)
call mpistop(
"cartesian thick EUV currently needs slab output")
567 call mpistop(
"thick EUV currently needs Cartesian dat_resolution output")
573 call mpistop(
"thick EUV currently needs x/y/z-aligned LOS")
580 double precision,
intent(in) :: emissivity,opacity,path_length
581 double precision,
intent(inout) :: intensity,tau
583 double precision :: dtau
585 if (path_length<=zero)
return
587 dtau=max(zero,opacity)*path_length
593 trim(emission_model)/=
'radio_ff' .and. &
598 logical,
intent(in) :: has_doppler,has_thick
601 if (has_doppler) num_outputs=num_outputs+1
602 if (has_thick .and. output_tau) num_outputs=num_outputs+1
603 if (has_thick .and. output_absorption_fraction) num_outputs=num_outputs+1
607 integer,
intent(in) :: nI1,nI2
608 double precision,
intent(in) :: EUV(nI1,nI2),unitv
609 double precision,
intent(inout) :: Dpl(nI1,nI2)
615 if (euv(ix1,ix2)/=zero)
then
616 dpl(ix1,ix2)=(dpl(ix1,ix2)/euv(ix1,ix2))*unitv
620 if (abs(dpl(ix1,ix2))<smalldouble) dpl(ix1,ix2)=zero
626 integer,
intent(in) :: nI1,nI2
627 double precision,
intent(in) :: EUV(nI1,nI2),EUVthin(nI1,nI2),smallflux
628 double precision,
intent(out) :: Absorption(nI1,nI2)
629 logical,
intent(in),
optional :: cap_to_one
632 logical :: cap_absorption
635 cap_absorption=.false.
636 if (
present(cap_to_one)) cap_absorption=cap_to_one
639 if (euvthin(ix1,ix2)>smallflux)
then
640 absorption(ix1,ix2)=max(zero,(euvthin(ix1,ix2)-euv(ix1,ix2))/euvthin(ix1,ix2))
641 if (cap_absorption) absorption(ix1,ix2)=min(one,absorption(ix1,ix2))
647 subroutine pack_euv_image_outputs(nI1,nI2,EUV,wI,smallflux,has_doppler,has_thick,Dpl,Tau,EUVthin,&
649 integer,
intent(in) :: nI1,nI2
650 double precision,
intent(in) :: EUV(nI1,nI2),smallflux
651 double precision,
intent(inout) :: wI(:,:,:)
652 logical,
intent(in) :: has_doppler,has_thick
653 double precision,
intent(in),
optional :: Dpl(nI1,nI2),Tau(nI1,nI2),EUVthin(nI1,nI2)
654 logical,
intent(in),
optional :: cap_absorption
657 double precision,
allocatable :: Absorption(:,:)
662 if (has_doppler)
then
663 if (.not.
present(dpl))
call mpistop(
"Doppler output requested without Doppler image")
667 if (has_thick .and. output_tau)
then
668 if (.not.
present(tau))
call mpistop(
"tau output requested without tau image")
672 if (has_thick .and. output_absorption_fraction)
then
673 if (.not.
present(euvthin))
call mpistop(
"absorption output requested without thin image")
674 allocate(absorption(ni1,ni2))
677 wi(:,:,iw)=absorption(:,:)
678 deallocate(absorption)
683 integer,
intent(out) :: pixel_batch_target,segment_batch_target,segment_comm_target
685 pixel_batch_target=max(1,radsyn_pixel_batch)
686 if (radsyn_segment_batch_factor>0)
then
687 segment_batch_target=max(128,radsyn_segment_batch_factor*pixel_batch_target)
689 segment_batch_target=max(128,int(min(dble(huge(segment_batch_target)),&
690 max(128.d0,radsyn_segment_memory_mb*1048576.d0/256.d0))))
692 segment_comm_target=max(128,radsyn_segment_comm_factor*pixel_batch_target)
696 double precision,
intent(in) :: tau
706 double precision,
intent(in) :: argument
708 if (argument<-700.d0)
then
710 else if (argument>700.d0)
then
718 double precision,
intent(in) :: exponent
720 if (exponent>300.d0)
then
722 else if (exponent<-300.d0)
then
730 double precision,
intent(in) :: temperature
731 integer,
intent(in) :: n_table
732 double precision,
intent(in) :: t_table(n_table),f_table(n_table)
733 logical,
intent(in) :: log_temperature,log_response
735 integer :: ilo,ihi,imid
736 double precision :: temp_lookup,response_lookup,flo,fhi
739 if (temperature<=zero)
return
740 if (log_temperature)
then
741 temp_lookup=log10(temperature)
743 temp_lookup=temperature
745 if (temp_lookup<t_table(1) .or. temp_lookup>t_table(n_table))
return
746 if (temp_lookup==t_table(n_table))
then
747 if (log_response)
then
748 response_lookup=log10(max(f_table(n_table),1.d-99))
750 response_lookup=f_table(n_table)
757 if (temp_lookup>=t_table(imid))
then
763 if (log_response)
then
764 flo=log10(max(f_table(ilo),1.d-99))
765 fhi=log10(max(f_table(ilo+1),1.d-99))
770 response_lookup=flo*(temp_lookup-t_table(ilo+1))/(t_table(ilo)-t_table(ilo+1))+&
771 fhi*(temp_lookup-t_table(ilo))/(t_table(ilo+1)-t_table(ilo))
774 if (log_response)
then
783 integer,
intent(in) :: ixI^L, ixO^L, n_table
784 double precision,
intent(in) :: Te(ixI^S),t_table(n_table),f_table(n_table)
785 double precision,
intent(inout) :: flux(ixI^S)
786 logical,
intent(in) :: log_temperature,log_response
789 double precision :: GT
791 {
do ix^db=ixomin^db,ixomax^db\}
793 flux(ix^d)=flux(ix^d)*gt
794 if (flux(ix^d)<zero) flux(ix^d)=zero
804 double precision,
intent(in) :: Te,Ne
805 double precision,
intent(out) :: x_HII,x_HeII,x_HeIII
807 double precision :: Pe,log_H21,log_He21,log_He32,log_He321
808 double precision :: logScaleHe,w_H21,term0,term1,term2,denHe
809 double precision,
parameter :: Xe_H21=13.6d0
810 double precision,
parameter :: Xe_He21=24.587d0
811 double precision,
parameter :: Xe_He32=54.416d0
821 log_h21=2.5d0*log10(te)-5040.d0*xe_h21/te-log10(pe)-0.48d0
822 log_he21=log10(4.d0)+2.5d0*log10(te)-5040.d0*xe_he21/te-log10(pe)-0.48d0
823 log_he32=2.5d0*log10(te)-5040.d0*xe_he32/te-log10(pe)-0.48d0
826 x_hii=w_h21/(1.d0+w_h21)
830 log_he321=log_he21+log_he32
831 logscalehe=max(
zero,log_he21,log_he321)
835 denhe=term0+term1+term2
848 double precision,
intent(in) :: nH,Te,rHe,Ne_guess
849 double precision,
intent(out) :: x_HII,x_HeII,x_HeIII
851 integer,
parameter :: max_iter=32
853 double precision :: Ne,Ne_lo,Ne_hi,Ne_new,residual,derivative
854 double precision :: e_He,de_HII_dNe,de_He_dNe
859 if (nh<=zero .or. te<=zero)
return
861 ne_lo=max(1.d-30*nh,1.d-100)
862 ne_hi=(1.d0+2.d0*rhe)*nh
863 ne=min(max(ne_guess,ne_lo),ne_hi)
867 e_he=x_heii+2.d0*x_heiii
868 residual=ne/nh-x_hii-rhe*e_he
869 if (abs(residual)<1.d-10)
exit
871 if (residual>zero)
then
879 de_hii_dne=-x_hii*(1.d0-x_hii)/ne
880 de_he_dne=(x_heii*(e_he-1.d0) &
881 +2.d0*x_heiii*(e_he-2.d0))/ne
882 derivative=1.d0/nh-de_hii_dne-rhe*de_he_dne
883 ne_new=ne-residual/derivative
884 if (.not.(ne_new>ne_lo .and. ne_new<ne_hi))
then
885 ne_new=0.5d0*(ne_lo+ne_hi)
899 integer,
intent(in) :: wl
900 integer,
intent(in) :: ixI^L, ixO^L
901 double precision,
intent(in) :: x(ixI^S,1:ndim)
902 double precision,
intent(in) :: w(ixI^S,1:nw)
904 double precision,
intent(out) :: kappa(ixI^S)
907 double precision :: pth(ixI^S),rho(ixI^S),Rfactor(ixI^S),Te(ixI^S)
908 double precision :: Ne(ixI^S),nH(ixI^S)
909 double precision :: wave_ratio,s_H1,s_He1,s_He2
910 double precision :: x_HII,x_HeII,x_HeIII,iz_H,iz_He,Rdummy
911 double precision :: N_H1,N_He1,N_He2
912 double precision,
parameter :: rHe_opacity=0.1d0
913 double precision,
parameter :: sigma_H1=5.16d-20, sigma_he1=9.25d-19, sigma_he2=7.17d-19
915 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
916 call fl%get_rho(w,x,ixi^l,ixo^l,rho)
917 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,rfactor)
919 {
do ix^db=ixomin^db,ixomax^db\}
920 if (rho(ix^d)>zero .and. rfactor(ix^d)>zero)
then
921 te(ix^d)=pth(ix^d)/(rho(ix^d)*rfactor(ix^d))*unit_temperature
928 call eos%get_ne_nH(ixi^l, ixo^l, w, x, ne, nh)
930 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
931 nh(ixo^s)=nh(ixo^s)*unit_numberdensity/1.d6
933 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
934 nh(ixo^s)=nh(ixo^s)*unit_numberdensity
937 wave_ratio=dble(wl)/171.d0
941 if (wl<=912) s_h1=wave_ratio**3*sigma_h1
942 if (wl<=504) s_he1=wave_ratio**2*sigma_he1
943 if (wl<=228) s_he2=wave_ratio**2.75d0*sigma_he2
946 {
do ix^db=ixomin^db,ixomax^db\}
947 if (te(ix^d)>zero .and. nh(ix^d)>zero)
then
948 select case (trim(eos%eos_type))
956 x_heii=iz_he*(1.d0-iz_he)
963 if (ne(ix^d)>zero)
then
965 x_hii,x_heii,x_heiii)
966 if (trim(eos%method)==
'analytic')
then
969 x_hii=min(one,max(zero,ne(ix^d)/nh(ix^d)))
971 x_hii=min(one,max(zero,ne(ix^d)/nh(ix^d) &
972 -eos%He_abundance*(x_heii+2.d0*x_heiii)))
976 rhe_opacity,ne(ix^d),x_hii,x_heii,x_heiii)
983 rhe_opacity,ne(ix^d),x_hii,x_heii,x_heiii)
986 n_h1=nh(ix^d)*(1.d0-x_hii)
987 n_he1=rhe_opacity*nh(ix^d)*(1.d0-x_heii-x_heiii)
988 n_he2=rhe_opacity*nh(ix^d)*x_heii
989 kappa(ix^d)=max(zero,n_h1*s_h1+n_he1*s_he1+n_he2*s_he2)
995 integer,
intent(in) :: igrid
996 integer,
intent(in) :: ixI^L, ixO^L
997 double precision,
intent(in) :: w(ixI^S,1:nw)
998 double precision,
intent(out) :: source(ixI^S)
1000 integer :: ix^D,idir,idirmin,idirmin0
1001 double precision :: current(ixI^S,7-2*ndir:3)
1003 if (.not.
allocated(iw_mag))
then
1004 call mpistop(
"emission_model=pseudo_current needs magnetic-field variables")
1009 call curlvector(w(ixi^s,iw_mag(1:ndir)),ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
1011 current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+ps(igrid)%J0(ixo^s,idirmin0:3)
1015 {
do ix^db=ixomin^db,ixomax^db\}
1017 source(ix^d)=source(ix^d)+current(ix^d,idir)**2
1025 integer,
intent(in) :: ixI^L, ixO^L
1026 double precision,
intent(in) :: x(ixI^S,1:ndim)
1027 double precision,
intent(in) :: w(ixI^S,1:nw)
1029 double precision,
intent(out) :: source(ixI^S),kappa(ixI^S)
1032 double precision :: pth(ixI^S),Te(ixI^S),Ne(ixI^S)
1033 double precision :: nH_dummy(ixI^S),gff
1035 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
1036 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
1037 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
1038 te(ixo^s)=pth(ixo^s)/(ne(ixo^s)*te(ixo^s))*unit_temperature
1039 call eos%get_ne_nH(ixi^l,ixo^l,w, x,ne,nh_dummy)
1041 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
1043 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
1048 {
do ix^db=ixomin^db,ixomax^db\}
1049 if (te(ix^d)>zero .and. ne(ix^d)>zero)
then
1050 if (te(ix^d)<2.d5)
then
1051 gff=18.2d0+1.5d0*log(te(ix^d))-log(radio_frequency)
1053 gff=24.5d0+log(te(ix^d))-log(radio_frequency)
1056 kappa(ix^d)=9.78d-3*ne(ix^d)**2*gff/(radio_frequency**2*te(ix^d)**1.5d0)
1057 source(ix^d)=te(ix^d)*kappa(ix^d)
1062 subroutine get_line_info(wl,ion,mass,logTe,line_center,spatial_px,spectral_px,sigma_PSF,width_slit)
1074 integer,
intent(in) :: wl
1075 integer,
intent(out) :: mass
1076 character(len=30),
intent(out) :: ion
1077 double precision,
intent(out) :: logTe,line_center,spatial_px,spectral_px
1078 double precision,
intent(out) :: sigma_PSF,width_slit
1148 line_center=1354.1d0
1150 spectral_px=12.98
d-3
1157 line_center=262.976d0
1166 line_center=263.765d0
1175 line_center=192.028d0
1184 line_center=255.113d0
1190 call mpistop(
"No information about this line")
1204 integer,
intent(in) :: wl
1205 integer,
intent(in) :: ixI^L, ixO^L
1206 double precision,
intent(in) :: x(ixI^S,1:ndim)
1207 double precision,
intent(in) :: w(ixI^S,1:nw)
1209 double precision,
intent(out) :: flux(ixI^S)
1212 double precision :: pth(ixI^S),rho(ixI^S),Rfactor(ixI^S),Te(ixI^S)
1213 double precision :: Ne(ixI^S),nH(ixI^S)
1215 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
1216 call fl%get_rho(w,x,ixi^l,ixo^l,rho)
1217 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,rfactor)
1223 call eos%get_ne_nH(ixi^l, ixo^l, w, x, ne, nh)
1231 flux(ixo^s)=ne(ixo^s)*nh(ixo^s)
1259 call mpistop(
"Unknown wavelength")
1271 integer,
intent(in) :: ixI^L,ixO^L
1272 integer,
intent(in) :: El,Eu
1273 double precision,
intent(in) :: x(ixI^S,1:ndim)
1274 double precision,
intent(in) :: w(ixI^S,nw)
1276 double precision,
intent(out) :: flux(ixI^S)
1278 integer :: ix^D,ixO^D
1280 double precision :: I0,kb,keV,dE,Ei
1281 double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
1282 double precision :: Ne(ixI^S),gff(ixI^S),fi(ixI^S)
1283 double precision :: EM(ixI^S)
1284 double precision :: nH_dummy(ixI^S)
1290 nume=floor((eu-el)/de)
1291 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
1292 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
1293 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
1296 call eos%get_ne_nH(ixi^l, ixo^l, w, x, ne, nh_dummy)
1299 em(ixo^s)=(ne(ixo^s))**2*1.d6
1302 em(ixo^s)=(ne(ixo^s))**2
1304 kbt(ixo^s)=kb*te(ixo^s)/kev
1309 {
do ix^db=ixomin^db,ixomax^db\}
1310 if (kbt(ix^d)>0.01*ei)
then
1311 if(kbt(ix^d)<ei) gff(ix^d)=(kbt(ix^d)/ei)**0.4
1312 fi(ix^d)=(em(ix^d)*gff(ix^d))*
exp_clamped(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
1317 flux(ixo^s)=flux(ixo^s)+fi(ixo^s)*de
1319 flux(ixo^s)=flux(ixo^s)*i0
1326 double precision,
intent(in) :: xbox^L
1328 double precision,
intent(out) :: eflux
1330 double precision :: dxb^D,xb^L
1331 integer :: iigrid,igrid,j
1332 integer :: ixO^L,ixI^L,ix^D
1333 double precision :: eflux_grid,eflux_pe
1335 ^d&iximin^d=
ixglo^d;
1336 ^d&iximax^d=
ixghi^d;
1337 ^d&ixomin^d=ixmlo^d;
1338 ^d&ixomax^d=ixmhi^d;
1340 do iigrid=1,igridstail; igrid=igrids(iigrid);
1344 call get_goes_flux_grid(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,ps(igrid)%dvolume(ixi^s),xbox^l,xb^l,fl,eflux_grid)
1345 eflux_pe=eflux_pe+eflux_grid
1347 call mpi_allreduce(eflux_pe,eflux,1,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
1354 integer,
intent(in) :: ixI^L,ixO^L
1355 double precision,
intent(in) :: x(ixI^S,1:ndim),dV(ixI^S)
1356 double precision,
intent(in) :: w(ixI^S,nw)
1357 double precision,
intent(in) :: xbox^L,xb^L
1359 double precision,
intent(out) :: eflux_grid
1361 integer :: ix^D,ixO^D,ixb^L
1362 integer :: iE,numE,j,inbox
1363 double precision :: I0,kb,keV,dE,Ei,El,Eu,A_cgs
1364 double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
1365 double precision :: Ne(ixI^S),EM(ixI^S)
1366 double precision :: nH_dummy(ixI^S)
1367 double precision :: gff,fi,erg_SI
1371 {
if (xbmin^d<xboxmax^d .and. xbmax^d>xboxmin^d) inbox=inbox+1\}
1373 if (inbox==ndim)
then
1375 ^d&ixbmin^d=ixomin^d;
1376 ^d&ixbmax^d=ixomax^d;
1377 {
if (xbmax^d>xboxmax^d) ixbmax^d=ixomax^d-ceiling((xbmax^d-xboxmax^d)/
dxlevel(^d))\}
1378 {
if (xbmin^d<xboxmin^d) ixbmin^d=ceiling((xboxmin^d-xbmin^d)/
dxlevel(^d))+ixomin^d\}
1385 el=const_h*const_c/(8.d0*a_cgs)/kev
1386 eu=const_h*const_c/(1.d0*a_cgs)/kev
1388 nume=floor((eu-el)/de)
1389 call fl%get_pthermal(w,x,ixi^l,ixb^l,pth)
1390 call fl%get_rho(w,x,ixi^l,ixb^l,ne)
1391 call fl%get_var_Rfactor(w,x,ixi^l,ixb^l,te)
1394 call eos%get_ne_nH(ixi^l, ixb^l, w, x, ne, nh_dummy)
1397 em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*(
unit_length*1.d2)**3
1400 em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*
unit_length**3
1402 kbt(ixb^s)=kb*te(ixb^s)/kev
1407 {
do ix^db=ixbmin^db,ixbmax^db\}
1408 if (kbt(ix^d)>1.d-2*ei)
then
1409 if(kbt(ix^d)<ei)
then
1410 gff=(kbt(ix^d)/ei)**0.4
1414 fi=(em(ix^d)*gff)*
exp_clamped(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
1415 eflux_grid=eflux_grid+fi*de*ei
1419 eflux_grid=eflux_grid*kev*erg_si
1428 integer,
intent(in) :: qunit
1430 character(20) :: datatype
1433 character (30) :: ion
1434 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1435 double precision :: xslit,arcsec
1437 datatype=
'spectrum_euv'
1442 if (
mype==0) print *,
'###################################################'
1445 if (
mype==0) print *,
'Systhesizing EUV spectrum (observed by IRIS).'
1446 case (263,264,192,255)
1447 if (
mype==0) print *,
'Systhesizing EUV spectrum (observed by Hinode/EIS).'
1449 call mpistop(
'Wrong wavelength!')
1453 call mpistop(
'Wrong spectrum window!')
1456 if (
mype==0)
write(*,
'(a,f8.3,a)')
' Wavelength: ',linecent,
' Angstrom'
1457 if (
mype==0) print *,
'Unit of EUV flux: DN s^-1 pixel^-1'
1461 write(*,
'(a,f5.3,a,f5.1,a)')
' Supposed pixel: ',wlrsl,
' Angstrom x ',spacersl*725.0,
' km'
1462 print *,
'Unit of wavelength: Angstrom (0.1 nm) '
1464 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
1466 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
1468 write(*,
'(a,f8.1,a)')
' Supposed width of slit: ',wslit*725.0,
' km'
1473 print *,
'Unit of wavelength: Angstrom (0.1 nm) '
1475 write(*,
'(a,f5.3,a,f5.1,a)')
' Pixel: ',wlrsl,
' Angstrom x ',spacersl*725.0,
' km'
1476 print *,
'Unit of length: arcsec (~725 km)'
1477 write(*,
'(a,f8.1,a)')
' Location of slit: xI1 = ',
location_slit,
' arcsec'
1478 write(*,
'(a,f8.1,a)')
' Width of slit: ',wslit,
' arcsec'
1481 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
1483 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
1485 write(*,
'(a,f8.1,a)')
' Location of slit: xI1 = ',
location_slit,
' Unit_length'
1486 write(*,
'(a,f8.1,a)')
' Width of slit: ',wslit*725.d0,
' km'
1489 if (
mype==0) print *,
'Direction of the slit: parallel to xI2 vector'
1493 call mpistop(
"EUV spectrum synthesis: support for sperical coordinates is to be added!")
1497 if (
mype==0) print *,
'###################################################'
1503 integer,
intent(in) :: qunit
1504 character(20),
intent(in) :: datatype
1507 integer :: numWL,numXS,iwL,ixS,numWI,numS
1508 double precision :: dwLg,xSmin,xSmax,wLmin,wLmax
1509 double precision,
allocatable :: wL(:),xS(:),dwL(:),dxS(:)
1510 double precision,
allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
1511 integer :: strtype,nstrb,nbb,nuni,nstr,bnx
1512 double precision :: qs,dxfirst,dxmid,lenstr
1514 integer :: iigrid,igrid,j,dir_loc
1515 double precision :: xbmin(1:ndim),xbmax(1:ndim)
1518 numwl=4*int((spectrum_window_max-spectrum_window_min)/(4.d0*dwlg))
1519 wlmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwlg*numwl/2
1520 wlmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwlg*numwl/2
1521 allocate(wl(numwl),dwl(numwl))
1524 wl(iwl)=wlmin+iwl*dwlg-half*dwlg
1527 select case(direction_slit)
1529 numxs=domain_nx1*2**(refine_max_level-1)
1534 strtype=stretch_type(1)
1535 nstrb=nstretchedblocks_baselevel(1)
1536 qs=qstretch_baselevel(1)
1537 if (mype==0) print *,
'Direction of the slit: x'
1539 numxs=domain_nx2*2**(refine_max_level-1)
1544 strtype=stretch_type(2)
1545 nstrb=nstretchedblocks_baselevel(2)
1546 qs=qstretch_baselevel(2)
1547 if (mype==0) print *,
'Direction of the slit: y'
1549 numxs=domain_nx3*2**(refine_max_level-1)
1554 strtype=stretch_type(3)
1555 nstrb=nstretchedblocks_baselevel(3)
1556 qs=qstretch_baselevel(3)
1557 if (mype==0) print *,
'Direction of the slit: z'
1559 call mpistop(
'Wrong direction_slit')
1562 allocate(xs(numxs),dxs(numxs),spectra(numwl,numxs),spectra_rc(numwl,numxs))
1564 allocate(wi(numwl,numxs,numwi))
1566 select case(strtype)
1568 dxs(:)=(xsmax-xsmin)/numxs
1570 xs(ixs)=xsmin+dxs(ixs)*(ixs-half)
1573 qs=qs**(one/2**(refine_max_level-1))
1574 dxfirst=(xsmax-xsmin)*(one-qs)/(one-qs**numxs)
1577 dxs(ixs)=dxfirst*qs**(ixs-1)
1578 xs(ixs)=dxs(1)/(one-qs)*(one-qs**(ixs-1))+half*dxs(ixs)
1584 lenstr=(xsmax-xsmin)/(2.d0+nuni*(one-qs)/(one-qs**nstr))
1585 dxfirst=(xsmax-xsmin)/(dble(nuni)+2.d0/(one-qs)*(one-qs**nstr))
1588 nstr=nstr*2**(refine_max_level-1)
1589 nuni=nuni*2**(refine_max_level-1)
1590 qs=qs**(one/2**(refine_max_level-1))
1591 dxfirst=lenstr*(one-qs)/(one-qs**nstr)
1592 dxmid=dxmid/2**(refine_max_level-1)
1594 if(nuni .gt. 0)
then
1595 do ixs=nstr+1,nstr+nuni
1597 xs(ixs)=lenstr+(dble(ixs)-0.5d0-nstr)*dxs(ixs)+xsmin
1602 dxs(ixs)=dxfirst*qs**(nstr-ixs)
1603 xs(ixs)=xsmin+lenstr-dxs(ixs)*half-dxfirst*(one-qs**(nstr-ixs))/(one-qs)
1606 do ixs=nstr+nuni+1,numxs
1607 dxs(ixs)=dxfirst*qs**(ixs-nstr-nuni-1)
1608 xs(ixs)=xsmax-lenstr+dxs(ixs)*half+dxfirst*(one-qs**(ixs-nstr-nuni-1))/(one-qs)
1611 call mpistop(
"unknown stretch type")
1614 if (los_phi==0 .and. los_theta==90 .and. direction_slit==2)
then
1617 else if (los_phi==0 .and. los_theta==90 .and. direction_slit==3)
then
1620 else if (los_phi==90 .and. los_theta==90 .and. direction_slit==1)
then
1623 else if (los_phi==90 .and. los_theta==90 .and. direction_slit==3)
then
1626 else if (los_theta==0 .and. direction_slit==1)
then
1629 else if (los_theta==0 .and. direction_slit==2)
then
1633 call mpistop(
'Wrong combination of LOS and slit direction!')
1636 if (dir_loc==1)
then
1637 if (location_slit>xprobmax1 .or. location_slit<xprobmin1)
then
1638 call mpistop(
'Wrong value for location_slit!')
1640 if(mype==0)
write(*,
'(a,f8.1,a)')
' Location of slit: x = ',location_slit,
' Unit_length'
1641 else if (dir_loc==2)
then
1642 if (location_slit>xprobmax2 .or. location_slit<xprobmin2)
then
1643 call mpistop(
'Wrong value for location_slit!')
1645 if(mype==0)
write(*,
'(a,f8.1,a)')
' Location of slit: y = ',location_slit,
' Unit_length'
1647 if (location_slit>xprobmax3 .or. location_slit<xprobmin3)
then
1648 call mpistop(
'Wrong value for location_slit!')
1650 if(mype==0)
write(*,
'(a,f8.1,a)')
' Location of slit: z = ',location_slit,
' Unit_length'
1655 do iigrid=1,igridstail; igrid=igrids(iigrid);
1656 ^d&xbmin(^d)=rnode(rpxmin^d_,igrid);
1657 ^d&xbmax(^d)=rnode(rpxmax^d_,igrid);
1658 if (location_slit>=xbmin(dir_loc) .and. location_slit<xbmax(dir_loc))
then
1664 call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
1665 mpi_sum,icomm,ierrmpi)
1668 if (spectra_rc(iwl,ixs)>smalldouble)
then
1669 wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
1676 call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
1678 deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
1685 integer,
intent(in) :: igrid,numWL,numXS,dir_loc
1687 double precision,
intent(in) :: wL(numWL),dwL(numWL)
1688 double precision,
intent(inout) :: spectra(numWL,numXS)
1690 integer :: direction_LOS
1691 integer :: ixO^L,ixI^L,ix^D,ixOnew
1692 double precision,
allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
1693 double precision :: wlc,wlwd
1696 double precision :: logTe,lineCent
1697 character (30) :: ion
1698 double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
1700 integer :: levelg,rft,ixSmin,ixSmax,iwL
1701 double precision :: flux_pix,dL
1703 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1705 if (los_phi==0 .and. los_theta==90)
then
1707 else if (los_phi==90 .and. los_theta==90)
then
1713 ^d&ixomin^d=ixmlo^d\
1714 ^d&ixomax^d=ixmhi^d\
1715 ^d&iximin^d=ixglo^d\
1716 ^d&iximax^d=ixghi^d\
1717 allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
1720 if (dir_loc==1)
then
1721 do ix1=ixomin1,ixomax1
1722 if (location_slit>=(ps(igrid)%x(ix^d,1)-
half*ps(igrid)%dx(ix^d,1)) .and. &
1723 location_slit<(ps(igrid)%x(ix^d,1)+
half*ps(igrid)%dx(ix^d,1)))
then
1729 else if (dir_loc==2)
then
1730 do ix2=ixomin2,ixomax2
1731 if (location_slit>=(ps(igrid)%x(ix^d,2)-
half*ps(igrid)%dx(ix^d,2)) .and. &
1732 location_slit<(ps(igrid)%x(ix^d,2)+
half*ps(igrid)%dx(ix^d,2)))
then
1739 do ix3=ixomin3,ixomax3
1740 if (location_slit>=(ps(igrid)%x(ix^d,3)-
half*ps(igrid)%dx(ix^d,3)) .and. &
1741 location_slit<(ps(igrid)%x(ix^d,3)+
half*ps(igrid)%dx(ix^d,3)))
then
1749 call get_euv(spectrum_wl,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
1750 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
1751 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
1752 v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
1753 call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
1754 call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
1755 te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
1758 levelg=ps(igrid)%level
1759 rft=2**(refine_max_level-levelg)
1761 {
do ix^d=ixomin^d,ixomax^d\}
1764 wlc=linecent*(1.d0+v(ix^d)*unit_velocity*1.d2/
const_c)
1766 wlc=linecent*(1.d0+v(ix^d)*unit_velocity/
const_c)
1768 wlwd=sqrt(
kb_cgs*te(ix^d)*unit_temperature/(mass*
mp_cgs))
1771 select case(direction_slit)
1773 ixsmin=(block_nx1*(node(pig1_,igrid)-1)+(ix1-ixomin1))*rft+1
1774 ixsmax=(block_nx1*(node(pig1_,igrid)-1)+(ix1-ixomin1+1))*rft
1776 ixsmin=(block_nx2*(node(pig2_,igrid)-1)+(ix2-ixomin2))*rft+1
1777 ixsmax=(block_nx2*(node(pig2_,igrid)-1)+(ix2-ixomin2+1))*rft
1779 ixsmin=(block_nx3*(node(pig3_,igrid)-1)+(ix3-ixomin3))*rft+1
1780 ixsmax=(block_nx3*(node(pig3_,igrid)-1)+(ix3-ixomin3+1))*rft
1783 select case(direction_los)
1785 dl=ps(igrid)%dx(ix^d,1)*unit_length
1787 dl=ps(igrid)%dx(ix^d,2)*unit_length
1789 dl=ps(igrid)%dx(ix^d,3)*unit_length
1791 if (si_unit) dl=dl*1.d2
1794 flux_pix=flux(ix^d)*wlrsl*dl*
exp_clamped(-(wl(iwl)-wlc)**2/(2*wlwd**2))/(sqrt(2*
dpi)*wlwd)
1796 flux_pix=flux_pix*wslit/spacersl
1797 spectra(iwl,ixsmin:ixsmax)=spectra(iwl,ixsmin:ixsmax)+flux_pix
1803 deallocate(flux,v,pth,te,rho)
1809 integer,
intent(in) :: qunit
1810 character(20),
intent(in) :: datatype
1813 integer :: numWL,numXS,iwL,ixS,numWI,ix^D
1814 double precision :: dwLg,dxSg,xSmin,xSmax,xScent,wLmin,wLmax
1815 double precision,
allocatable :: wL(:),xS(:),dwL(:),dxS(:)
1816 double precision,
allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
1817 double precision :: vec_cor(1:3),xI_cor(1:2)
1818 double precision :: res,r_loc,r_max
1821 character (30) :: ion
1822 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1823 double precision :: unitv,arcsec,RHESSI_rsl,pixel
1824 integer :: iigrid,igrid,i,j,numS
1825 double precision :: xLmin,xLmax,xslit
1836 xsmin=-abs(xprobmax1)
1837 xsmax=abs(xprobmax1)
1840 if (ix1==1) vec_cor(1)=xprobmin1
1841 if (ix1==2) vec_cor(1)=xprobmax1
1843 if (ix2==1) vec_cor(2)=xprobmin2
1844 if (ix2==2) vec_cor(2)=xprobmax2
1846 if (ix3==1) vec_cor(3)=xprobmin3
1847 if (ix3==2) vec_cor(3)=xprobmax3
1849 r_loc=(vec_cor(1)-x_origin(1))**2
1850 r_loc=r_loc+(vec_cor(2)-x_origin(2))**2
1851 r_loc=r_loc+(vec_cor(3)-x_origin(3))**2
1853 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1856 r_max=max(r_max,r_loc)
1860 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1864 xsmin=min(xsmin,xi_cor(2))
1865 xsmax=max(xsmax,xi_cor(2))
1876 xscent=(xsmin+xsmax)/2.d0
1880 arcsec=7.25d5/unit_length
1882 arcsec=7.25d7/unit_length
1884 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1885 dxsg=spacersl*arcsec
1886 numxs=ceiling((xsmax-xscent)/dxsg)
1887 xsmin=xscent-numxs*dxsg
1888 xsmax=xscent+numxs*dxsg
1891 numwl=2*int((spectrum_window_max-spectrum_window_min)/(2.d0*dwlg))
1892 wlmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwlg*numwl/2
1893 wlmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwlg*numwl/2
1894 allocate(wl(numwl),dwl(numwl),xs(numxs),dxs(numxs))
1896 allocate(wi(numwl,numxs,numwi),spectra(numwl,numxs),spectra_rc(numwl,numxs))
1898 wl(iwl)=wlmin+iwl*dwlg-half*dwlg
1902 xs(ixs)=xsmin+dxsg*(ixs-half)
1908 do iigrid=1,igridstail; igrid=igrids(iigrid);
1910 if (ix1==1) vec_cor(1)=rnode(rpxmin1_,igrid)
1911 if (ix1==2) vec_cor(1)=rnode(rpxmax1_,igrid)
1913 if (ix2==1) vec_cor(2)=rnode(rpxmin2_,igrid)
1914 if (ix2==2) vec_cor(2)=rnode(rpxmax2_,igrid)
1916 if (ix3==1) vec_cor(3)=rnode(rpxmin3_,igrid)
1917 if (ix3==2) vec_cor(3)=rnode(rpxmax3_,igrid)
1919 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1923 xlmin=min(xlmin,xi_cor(1))
1924 xlmax=max(xlmax,xi_cor(1))
1930 if (activate_unit_arcsec)
then
1931 xslit=location_slit*arcsec
1935 if (xslit>=xlmin-wslit*arcsec .and. xslit<=xlmax+wslit*arcsec)
then
1941 call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
1942 mpi_sum,icomm,ierrmpi)
1945 if (spectra_rc(iwl,ixs)>smalldouble)
then
1946 wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
1953 if (activate_unit_arcsec)
then
1958 call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
1960 deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
1966 integer,
intent(in) :: igrid,numWL,numXS
1967 double precision,
intent(in) :: wL(numWL),xS(numXS)
1968 double precision,
intent(in) :: dwLg,dxSg
1969 double precision,
intent(inout) :: spectra(numWL,numXS)
1972 integer :: ixO^L,ixI^L,ix^D,ixOnew,j
1973 double precision,
allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
1974 double precision :: wlc,wlwd,res,dst_slit,xslit,arcsec
1975 double precision :: vloc(1:3),xloc(1:3),dxloc(1:3),xIloc(1:2),dxIloc(1:2)
1976 integer :: nSubC^D,iSubC^D,iwL,ixS,ixSmin,ixSmax,iwLmin,iwLmax,nwL
1977 double precision :: slit_width,dxSubC^D,xerf^L,fluxSubC
1978 double precision :: xSubC(1:3),xCent(1:2)
1981 double precision :: logTe,lineCent
1982 character (30) :: ion
1983 double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
1984 double precision :: sigma_wl,sigma_xs,factor
1987 arcsec=7.25d5/unit_length
1989 arcsec=7.25d7/unit_length
1991 if (activate_unit_arcsec)
then
1992 xslit=location_slit*arcsec
1997 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1999 ^d&ixomin^d=ixmlo^d\
2000 ^d&ixomax^d=ixmhi^d\
2001 ^d&iximin^d=ixglo^d\
2002 ^d&iximax^d=ixghi^d\
2003 allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
2005 call get_euv(spectrum_wl,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
2006 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
2007 call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
2008 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
2009 call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
2010 te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
2011 {
do ix^d=ixomin^d,ixomax^d\}
2013 vloc(j)=ps(igrid)%w(ix^d,iw_mom(j))/rho(ix^d)
2021 slit_width=wslit*arcsec
2022 sigma_wl=sigma_psf*dwlg
2023 sigma_xs=sigma_psf*dxsg
2024 {
do ix^d=ixomin^d,ixomax^d\}
2025 if (flux(ix^d)>smalldouble)
then
2026 xloc(1:3)=ps(igrid)%x(ix^d,1:3)
2027 dxloc(1:3)=ps(igrid)%dx(ix^d,1:3)
2031 if (xiloc(1)>=xslit-half*(slit_width+dxiloc(1)) .and. &
2032 xiloc(1)<=xslit+half*(slit_width+dxiloc(1)))
then
2034 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi1(^d))/(slit_width/16.d0)));
2035 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi2(^d))/(dxsg/4.d0)));
2036 ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
2039 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length*1.d2/dxsg/dxsg
2040 wlc=linecent*(1.d0+v(ix^d)*unit_velocity*1.d2/const_c)
2042 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length/dxsg/dxsg
2043 wlc=linecent*(1.d0+v(ix^d)*unit_velocity/const_c)
2045 wlwd=sqrt(kb_cgs*te(ix^d)*unit_temperature/(mass*mp_cgs))
2046 wlwd=wlwd*linecent/const_c
2048 {
do isubc^d=1,nsubc^d\}
2049 ^d&xsubc(^d)=xloc(^d)-half*dxloc(^d)+(isubc^d-half)*dxsubc^d;
2051 dst_slit=abs(xcent(1)-xslit)
2052 if (dst_slit<=half*slit_width)
then
2053 ixs=floor((xcent(2)-(xs(1)-half*dxsg))/dxsg)+1
2055 ixsmax=min(ixs+3,numxs)
2056 iwl=floor((wlc-(wl(1)-half*dwlg))/dwlg)+1
2057 nwl=3*ceiling(wlwd/dwlg+1)
2058 iwlmin=max(1,iwl-nwl)
2059 iwlmax=min(iwl+nwl,numwl)
2061 do iwl=iwlmin,iwlmax
2062 do ixs=ixsmin,ixsmax
2063 xerfmin1=(wl(iwl)-half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
2064 xerfmax1=(wl(iwl)+half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
2065 xerfmin2=(xs(ixs)-half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
2066 xerfmax2=(xs(ixs)+half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
2067 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
2068 spectra(iwl,ixs)=spectra(iwl,ixs)+fluxsubc*factor
2078 deallocate(flux,v,pth,te)
2086 integer,
intent(in) :: qunit
2088 character(20) :: datatype
2091 character (30) :: ion
2092 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
2093 double precision :: t0,t1
2096 datatype=
'image_euv'
2101 print *,
'###################################################'
2102 print *,
'Systhesizing EUV image'
2103 write(*,
'(a,f8.3,a)')
' Wavelength: ',linecent,
' Angstrom'
2104 print *,
'Unit of EUV flux: DN s^-1 pixel^-1'
2109 call mpistop(
'EUV dat-resolution needs Cartesian or spherical native rays')
2111 print *,
'Data-resolution image requested; native output pixel sizes are reported below.'
2113 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2115 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2129 call mpistop(
'ERROR: Wrong LOS for synthesizing emission!')
2133 write(*,
'(a,f7.1,a,f7.1,a,f5.1,a,f5.1,a)')
' Pixel: ',spacersl*725.0,
' km x ',spacersl*725.0,
' km (', &
2134 spacersl,
' arcsec x ', spacersl,
' arcsec)'
2136 print *,
'Unit of length: arcsec (~725 km)'
2139 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2141 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2147 '] of the simulation box is located at [X=0,Y=0] of the image'
2150 if (
mype==0)
write(*,
'(a,f6.3,f8.3,f8.3,a)')
' Mapping: R=0 of the simulation box is located at [X=0,Y=0] of the image'
2153 call mpistop(
"EUV synthesis: this coordinate is not supported!")
2158 if (
mype==0) print *,
'time comsuming: ',t1-t0,
' s'
2159 if (
mype==0) print *,
'###################################################'
2166 integer,
intent(in) :: qunit
2168 character(20) :: datatype
2169 double precision :: RHESSI_rsl
2170 double precision :: t0,t1
2173 datatype=
'image_sxr'
2178 print *,
'###################################################'
2179 print *,
'Systhesizing SXR image (observed at 1 AU).'
2186 print *,
'Unit of SXR flux: photons cm^-2 s^-1 pixel^-1'
2187 write(*,
'(a,f5.1,a,f5.1,a,f5.1,a,f5.1,a)')
' Supposed Pixel: ',rhessi_rsl*0.725,
' Mm x ',rhessi_rsl*0.725, &
2188 ' Mm (', rhessi_rsl,
' arcsec x ', rhessi_rsl,
' arcsec)'
2190 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2192 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2202 call mpistop(
'ERROR: Wrong LOS for synthesizing emission!')
2206 print *,
'Unit of SXR flux: photons cm^-2 s^-1 pixel^-1'
2207 write(*,
'(a,f5.1,a,f5.1,a,f5.1,a,f5.1,a)')
' Pixel: ',rhessi_rsl*0.725,
' Mm x ',rhessi_rsl*0.725, &
2208 ' Mm (', rhessi_rsl,
' arcsec x ', rhessi_rsl,
' arcsec)'
2210 print *,
'Unit of length: arcsec (~725 km)'
2213 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2215 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2221 '] of the simulation box is located at [X=0,Y=0] of the image'
2224 if (
mype==0)
write(*,
'(a,f6.3,f8.3,f8.3,a)')
' Mapping: R=0 of the simulation box is located at [X=0,Y=0] of the image'
2227 call mpistop(
"SXR synthesis: this coordinate is not supported!")
2232 if (
mype==0) print *,
'time comsuming:',t1-t0
2233 if (
mype==0) print *,
'###################################################'
2240 integer,
intent(in) :: qunit
2242 character(20) :: datatype
2243 double precision :: LASCO_rsl
2245 if (
mype==0) print *,
'###################################################'
2249 if (
mype==0) print *,
'Systhesizing white light image (observed by LASCO/C1).'
2252 if (
mype==0) print *,
'Systhesizing white light image (observed by LASCO/C2).'
2255 if (
mype==0) print *,
'Systhesizing white light image (observed by LASCO/C3).'
2257 call mpistop(
'Whitelight synthesis: instrument is not supported!')
2260 if (
mype==0)
write(*,
'(a,f5.1,a,f5.1,a,f5.1,a,f5.1,a)')
' Pixel: ',lasco_rsl*0.725,
' Mm x ',lasco_rsl*0.725,
' Mm (', &
2261 lasco_rsl,
' arcsec x ', lasco_rsl,
' arcsec) '
2262 if (
mype==0) print *,
'Unit of white light flux: average Sun brightness'
2264 datatype=
'image_whitelight'
2269 print *,
'Unit of length: arcsec (~725 km)'
2272 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2274 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2280 if (
mype==0)
write(*,
'(a,f6.3,f8.3,f8.3,a)')
' Mapping: R=0 of the simulation box is located at [X=0,Y=0] of the image'
2283 call mpistop(
"Whitelight synthesis: this coordinate is not supported!")
2286 if (
mype==0) print *,
'###################################################'
2291 EUV,Dpl,nOut1,nOut2,xOut1,xOut2,&
2292 dxOut1,dxOut2,wOut,numWOut,Tau,EUVthin)
2296 integer,
intent(in) :: nSrc1,nSrc2
2297 double precision,
intent(in) :: xSrc1(nSrc1),xSrc2(nSrc2)
2298 double precision,
intent(in) :: dxSrc1(nSrc1),dxSrc2(nSrc2)
2299 double precision,
intent(in) :: EUV(nSrc1,nSrc2),Dpl(nSrc1,nSrc2)
2300 integer,
intent(out) :: nOut1,nOut2,numWOut
2301 double precision,
allocatable,
intent(out) :: xOut1(:),xOut2(:),dxOut1(:),dxOut2(:)
2302 double precision,
allocatable,
intent(out) :: wOut(:,:,:)
2303 double precision,
intent(in),
optional :: Tau(nSrc1,nSrc2),EUVthin(nSrc1,nSrc2)
2305 integer :: mass,ixS1,ixS2,ixP1,ixP2,ixC1,ixC2,iw
2306 integer :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
2307 character(30) :: ion
2308 double precision :: logTe,lineCent,spaceRsl,wlRsl,sigma_PSF,wslit
2309 double precision :: arcsec,dxInst,xMin1,xMax1,xMin2,xMax2,xCent1,xCent2
2310 double precision :: sigma0,xerfmin1,xerfmax1,xerfmin2,xerfmax2
2311 double precision :: factor,weightSum,weightNorm,thinVal,tauVal
2312 double precision,
allocatable :: dplNum(:,:),thinOut(:,:),tauOut(:,:),tauWeight(:,:)
2320 dxinst=spacersl*arcsec
2321 if (dxinst<=
zero)
call mpistop(
"instrument_postprocess has non-positive pixel size")
2323 xmin1=minval(xsrc1-
half*dxsrc1)
2324 xmax1=maxval(xsrc1+
half*dxsrc1)
2325 xmin2=minval(xsrc2-
half*dxsrc2)
2326 xmax2=maxval(xsrc2+
half*dxsrc2)
2327 xcent1=
half*(xmin1+xmax1)
2328 xcent2=
half*(xmin2+xmax2)
2329 nout1=16*max(1,ceiling((xmax1-xmin1)/(16.d0*dxinst)))
2330 nout2=16*max(1,ceiling((xmax2-xmin2)/(16.d0*dxinst)))
2331 xmin1=xcent1-
half*dble(nout1)*dxinst
2332 xmin2=xcent2-
half*dble(nout2)*dxinst
2334 allocate(xout1(nout1),xout2(nout2),dxout1(nout1),dxout2(nout2))
2336 xout1(ixp1)=xmin1+dxinst*(dble(ixp1)-
half)
2340 xout2(ixp2)=xmin2+dxinst*(dble(ixp2)-
half)
2345 if (
present(tau) .and.
output_tau) numwout=numwout+1
2347 allocate(wout(nout1,nout2,numwout),dplnum(nout1,nout2))
2350 if (
present(euvthin))
then
2351 allocate(thinout(nout1,nout2))
2355 allocate(tauout(nout1,nout2),tauweight(nout1,nout2))
2360 sigma0=sigma_psf*dxinst
2365 if (
present(euvthin)) thinval=euvthin(ixs1,ixs2)
2366 if (
present(tau)) tauval=tau(ixs1,ixs2)
2370 ixc1=floor((xsrc1(ixs1)-(xout1(1)-
half*dxinst))/dxinst)+1
2371 ixc2=floor((xsrc2(ixs2)-(xout2(1)-
half*dxinst))/dxinst)+1
2372 ixpmin1=max(1,ixc1-3)
2373 ixpmax1=min(nout1,ixc1+3)
2374 ixpmin2=max(1,ixc2-3)
2375 ixpmax2=min(nout2,ixc2+3)
2378 do ixp1=ixpmin1,ixpmax1
2379 do ixp2=ixpmin2,ixpmax2
2380 xerfmin1=((xout1(ixp1)-
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2381 xerfmax1=((xout1(ixp1)+
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2382 xerfmin2=((xout2(ixp2)-
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2383 xerfmax2=((xout2(ixp2)+
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2384 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
2385 weightsum=weightsum+factor
2388 if (weightsum<=
zero) cycle
2390 do ixp1=ixpmin1,ixpmax1
2391 do ixp2=ixpmin2,ixpmax2
2392 xerfmin1=((xout1(ixp1)-
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2393 xerfmax1=((xout1(ixp1)+
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2394 xerfmin2=((xout2(ixp2)-
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2395 xerfmax2=((xout2(ixp2)+
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2396 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
2397 weightnorm=factor/weightsum
2398 wout(ixp1,ixp2,1)=wout(ixp1,ixp2,1)+euv(ixs1,ixs2)*weightnorm
2399 dplnum(ixp1,ixp2)=dplnum(ixp1,ixp2)+euv(ixs1,ixs2)*dpl(ixs1,ixs2)*weightnorm
2400 if (
present(euvthin)) thinout(ixp1,ixp2)=thinout(ixp1,ixp2)+thinval*weightnorm
2402 tauout(ixp1,ixp2)=tauout(ixp1,ixp2)+tauval*weightnorm
2403 tauweight(ixp1,ixp2)=tauweight(ixp1,ixp2)+weightnorm
2413 wout(ixp1,ixp2,2)=dplnum(ixp1,ixp2)/wout(ixp1,ixp2,1)
2415 wout(ixp1,ixp2,2)=
zero
2424 if (tauweight(ixp1,ixp2)>
zero)
then
2425 wout(ixp1,ixp2,iw)=tauout(ixp1,ixp2)/tauweight(ixp1,ixp2)
2427 wout(ixp1,ixp2,iw)=
zero
2437 wout(ixp1,ixp2,iw)=min(
one,max(
zero,(thinout(ixp1,ixp2)-wout(ixp1,ixp2,1))/thinout(ixp1,ixp2)))
2439 wout(ixp1,ixp2,iw)=
zero
2446 write(*,
'(a,2(i8,1x),a,2(i8,1x),a,1pe12.5)') &
2447 ' instrument_postprocess EUV grid src/out: ',nsrc1,nsrc2,
' -> ',nout1,nout2,
' dx=',dxinst
2451 if (
allocated(thinout))
deallocate(thinout)
2452 if (
allocated(tauout))
deallocate(tauout,tauweight)
2456 Bright,nOut1,nOut2,xOut1,xOut2,&
2457 dxOut1,dxOut2,wOut,numWOut,Tau,BrightThin)
2460 integer,
intent(in) :: nSrc1,nSrc2
2461 double precision,
intent(in) :: xSrc1(nSrc1),xSrc2(nSrc2)
2462 double precision,
intent(in) :: dxSrc1(nSrc1),dxSrc2(nSrc2)
2463 double precision,
intent(in) :: Bright(nSrc1,nSrc2)
2464 integer,
intent(out) :: nOut1,nOut2,numWOut
2465 double precision,
allocatable,
intent(out) :: xOut1(:),xOut2(:),dxOut1(:),dxOut2(:)
2466 double precision,
allocatable,
intent(out) :: wOut(:,:,:)
2467 double precision,
intent(in),
optional :: Tau(nSrc1,nSrc2),BrightThin(nSrc1,nSrc2)
2469 integer :: ixS1,ixS2,ixP1,ixP2,ixC1,ixC2,iw,nStencil
2470 integer :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
2471 double precision :: arcsec,beamPixel,beamSigma,xMin1,xMax1,xMin2,xMax2,xCent1,xCent2
2472 double precision :: distance1,distance2,weight,cellArea,thinVal,tauVal
2473 double precision,
allocatable :: norm(:,:),thinOut(:,:),tauOut(:,:),tauNorm(:,:)
2486 if (beamsigma<=zero .or. beampixel<=zero)
then
2487 call mpistop(
"radio beam postprocess has non-positive beam or pixel size")
2490 xmin1=minval(xsrc1-half*dxsrc1)
2491 xmax1=maxval(xsrc1+half*dxsrc1)
2492 xmin2=minval(xsrc2-half*dxsrc2)
2493 xmax2=maxval(xsrc2+half*dxsrc2)
2494 xcent1=half*(xmin1+xmax1)
2495 xcent2=half*(xmin2+xmax2)
2496 nout1=16*max(1,ceiling((xmax1-xmin1)/(16.d0*beampixel)))
2497 nout2=16*max(1,ceiling((xmax2-xmin2)/(16.d0*beampixel)))
2498 xmin1=xcent1-half*dble(nout1)*beampixel
2499 xmin2=xcent2-half*dble(nout2)*beampixel
2501 allocate(xout1(nout1),xout2(nout2),dxout1(nout1),dxout2(nout2))
2503 xout1(ixp1)=xmin1+beampixel*(dble(ixp1)-half)
2504 dxout1(ixp1)=beampixel
2507 xout2(ixp2)=xmin2+beampixel*(dble(ixp2)-half)
2508 dxout2(ixp2)=beampixel
2512 if (
present(tau) .and.
output_tau) numwout=numwout+1
2514 allocate(wout(nout1,nout2,numwout),norm(nout1,nout2))
2517 if (
present(brightthin))
then
2518 allocate(thinout(nout1,nout2))
2522 allocate(tauout(nout1,nout2),taunorm(nout1,nout2))
2527 nstencil=max(3,ceiling(4.d0*beamsigma/beampixel)+1)
2532 if (
present(brightthin)) thinval=brightthin(ixs1,ixs2)
2533 if (
present(tau)) tauval=tau(ixs1,ixs2)
2534 if (abs(bright(ixs1,ixs2))<=smalldouble .and. abs(thinval)<=smalldouble .and. &
2535 abs(tauval)<=smalldouble) cycle
2537 ixc1=floor((xsrc1(ixs1)-(xout1(1)-half*beampixel))/beampixel)+1
2538 ixc2=floor((xsrc2(ixs2)-(xout2(1)-half*beampixel))/beampixel)+1
2539 ixpmin1=max(1,ixc1-nstencil)
2540 ixpmax1=min(nout1,ixc1+nstencil)
2541 ixpmin2=max(1,ixc2-nstencil)
2542 ixpmax2=min(nout2,ixc2+nstencil)
2543 cellarea=max(smalldouble,dxsrc1(ixs1)*dxsrc2(ixs2))
2545 do ixp1=ixpmin1,ixpmax1
2546 distance1=xout1(ixp1)-xsrc1(ixs1)
2547 do ixp2=ixpmin2,ixpmax2
2548 distance2=xout2(ixp2)-xsrc2(ixs2)
2549 weight=
exp_clamped(-half*(distance1**2+distance2**2)/beamsigma**2)*cellarea
2550 wout(ixp1,ixp2,1)=wout(ixp1,ixp2,1)+bright(ixs1,ixs2)*weight
2551 norm(ixp1,ixp2)=norm(ixp1,ixp2)+weight
2552 if (
present(brightthin)) thinout(ixp1,ixp2)=thinout(ixp1,ixp2)+thinval*weight
2554 tauout(ixp1,ixp2)=tauout(ixp1,ixp2)+tauval*weight
2555 taunorm(ixp1,ixp2)=taunorm(ixp1,ixp2)+weight
2564 if (norm(ixp1,ixp2)>zero)
then
2565 wout(ixp1,ixp2,1)=wout(ixp1,ixp2,1)/norm(ixp1,ixp2)
2566 if (
present(brightthin)) thinout(ixp1,ixp2)=thinout(ixp1,ixp2)/norm(ixp1,ixp2)
2568 wout(ixp1,ixp2,1)=zero
2569 if (
present(brightthin)) thinout(ixp1,ixp2)=zero
2579 if (taunorm(ixp1,ixp2)>zero)
then
2580 wout(ixp1,ixp2,iw)=tauout(ixp1,ixp2)/taunorm(ixp1,ixp2)
2582 wout(ixp1,ixp2,iw)=zero
2591 if (thinout(ixp1,ixp2)>smalldouble)
then
2592 wout(ixp1,ixp2,iw)=min(one,max(zero,(thinout(ixp1,ixp2)-wout(ixp1,ixp2,1))/thinout(ixp1,ixp2)))
2594 wout(ixp1,ixp2,iw)=zero
2601 write(*,
'(a,2(i8,1x),a,2(i8,1x),a,2(1pe12.5,1x))') &
2602 ' radio_beam_postprocess grid src/out: ',nsrc1,nsrc2,
' -> ',nout1,nout2,&
2607 if (
allocated(thinout))
deallocate(thinout)
2608 if (
allocated(tauout))
deallocate(tauout,taunorm)
2617 integer,
intent(in) :: qunit
2618 character(20),
intent(in) :: datatype
2621 double precision :: dx^D
2622 integer :: numX^D,ix^D
2623 double precision,
allocatable :: EUV(:,:),EUVs(:,:),Dpl(:,:),Dpls(:,:)
2624 double precision,
allocatable :: EUVthin(:,:),Tau(:,:)
2625 double precision,
allocatable :: SXR(:,:),SXRs(:,:),wI(:,:,:)
2626 double precision,
allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:),dxIi
2627 integer :: numXI1,numXI2,numSI,numWI,iw
2628 double precision :: xI^L
2629 integer :: iigrid,igrid,i,j
2630 double precision,
allocatable :: xIF1(:),xIF2(:),dxIF1(:),dxIF2(:)
2631 double precision,
allocatable :: xIP1(:),xIP2(:),dxIP1(:),dxIP2(:),wIP(:,:,:)
2632 integer :: nXIF1,nXIF2
2633 integer :: nXIP1,nXIP2,numWIP
2634 double precision :: xIF^L
2635 double precision :: vec_cor(1:3),xI_cor(1:2),dxDDA,xIcent1,xIcent2
2637 double precision :: unitv,arcsec,RHESSI_rsl,length_to_km
2638 integer :: strtype^D,nstrb^D,nbb^D,nuni^D,nstr^D,bnx^D
2639 double precision :: qs^D,dxfirst^D,dxmid^D,lenstr^D
2640 logical :: has_doppler_output,has_thick_output
2653 xicent1=
half*(xifmin1+xifmax1)
2654 xicent2=
half*(xifmin2+xifmax2)
2655 nxif1=max(1,ceiling((xifmax1-xifmin1)/dxdda))
2656 nxif2=max(1,ceiling((xifmax2-xifmin2)/dxdda))
2657 xifmin1=xicent1-
half*dble(nxif1)*dxdda
2658 xifmax1=xicent1+
half*dble(nxif1)*dxdda
2659 xifmin2=xicent2-
half*dble(nxif2)*dxdda
2660 xifmax2=xicent2+
half*dble(nxif2)*dxdda
2671 if (
mype==0)
write(*,
'(a,a,a,1pe12.5,a,2(i8,1x))') &
2673 ' image-plane dx=',dxdda,
' n=',nxif1,nxif2
2676 if (ix1==1) vec_cor(1)=xprobmin1
2677 if (ix1==2) vec_cor(1)=xprobmax1
2679 if (ix2==1) vec_cor(2)=xprobmin2
2680 if (ix2==2) vec_cor(2)=xprobmax2
2682 if (ix3==1) vec_cor(3)=xprobmin3
2683 if (ix3==2) vec_cor(3)=xprobmax3
2685 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
2691 xifmin1=min(xifmin1,xi_cor(1))
2692 xifmax1=max(xifmax1,xi_cor(1))
2693 xifmin2=min(xifmin2,xi_cor(2))
2694 xifmax2=max(xifmax2,xi_cor(2))
2700 xicent1=
half*(xifmin1+xifmax1)
2701 xicent2=
half*(xifmin2+xifmax2)
2702 nxif1=max(1,ceiling((xifmax1-xifmin1)/dxdda))
2703 nxif2=max(1,ceiling((xifmax2-xifmin2)/dxdda))
2704 xifmin1=xicent1-
half*dble(nxif1)*dxdda
2705 xifmax1=xicent1+
half*dble(nxif1)*dxdda
2706 xifmin2=xicent2-
half*dble(nxif2)*dxdda
2707 xifmax2=xicent2+
half*dble(nxif2)*dxdda
2718 if (
mype==0)
write(*,
'(a,a,a,1pe12.5,a,2(i8,1x))') &
2720 ' image-plane dx=',dxdda,
' n=',nxif1,nxif2
2738 if (
mype==0)
write(*,
'(a)')
' LOS vector: [-1.00 0.00 0.00]'
2739 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [ 0.00 1.00 0.00]'
2740 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 0.00 1.00]'
2758 if (
mype==0)
write(*,
'(a)')
' LOS vector: [ 0.00 -1.00 0.00]'
2759 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [-1.00 0.00 0.00]'
2760 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 0.00 1.00]'
2778 if (
mype==0)
write(*,
'(a)')
' LOS vector: [ 0.00 0.00 -1.00]'
2779 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [ 1.00 0.00 0.00]'
2780 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 1.00 0.00]'
2782 allocate(xif1(nxif1),xif2(nxif2),dxif1(nxif1),dxif2(nxif2))
2785 select case(strtype1)
2787 dxif1(:)=(xifmax1-xifmin1)/nxif1
2789 xif1(ix1)=xifmin1+dxif1(ix1)*(ix1-
half)
2793 dxfirst1=(xifmax1-xifmin1)*(
one-qs1)/(
one-qs1**nxif1)
2796 dxif1(ix1)=dxfirst1*qs1**(ix1-1)
2797 xif1(ix1)=dxif1(1)/(
one-qs1)*(
one-qs1**(ix1-1))+
half*dxif1(ix1)
2802 nuni1=nbb1-nstrb1*bnx1
2803 lenstr1=(xifmax1-xifmin1)/(2.d0+nuni1*(
one-qs1)/(
one-qs1**nstr1))
2804 dxfirst1=(xifmax1-xifmin1)/(dble(nuni1)+2.d0/(
one-qs1)*(
one-qs1**nstr1))
2810 dxfirst1=lenstr1*(
one-qs1)/(
one-qs1**nstr1)
2813 if(nuni1 .gt. 0)
then
2814 do ix1=nstr1+1,nstr1+nuni1
2816 xif1(ix1)=lenstr1+(dble(ix1)-0.5d0-nstr1)*dxif1(ix1)+xifmin1
2821 dxif1(ix1)=dxfirst1*qs1**(nstr1-ix1)
2822 xif1(ix1)=xifmin1+lenstr1-dxif1(ix1)*
half-dxfirst1*(
one-qs1**(nstr1-ix1))/(
one-qs1)
2825 do ix1=nstr1+nuni1+1,nxif1
2826 dxif1(ix1)=dxfirst1*qs1**(ix1-nstr1-nuni1-1)
2827 xif1(ix1)=xifmax1-lenstr1+dxif1(ix1)*
half+dxfirst1*(
one-qs1**(ix1-nstr1-nuni1-1))/(
one-qs1)
2830 call mpistop(
"unknown stretch type")
2833 select case(strtype2)
2835 dxif2(:)=(xifmax2-xifmin2)/nxif2
2837 xif2(ix2)=xifmin2+dxif2(ix2)*(ix2-
half)
2841 dxfirst2=(xifmax2-xifmin2)*(
one-qs2)/(
one-qs2**nxif2)
2844 dxif2(ix2)=dxfirst2*qs2**(ix2-1)
2845 xif2(ix2)=dxif2(1)/(
one-qs2)*(
one-qs2**(ix2-1))+
half*dxif2(ix2)
2850 nuni2=nbb2-nstrb2*bnx2
2851 lenstr2=(xifmax2-xifmin2)/(2.d0+nuni2*(
one-qs2)/(
one-qs2**nstr2))
2852 dxfirst2=(xifmax2-xifmin2)/(dble(nuni2)+2.d0/(
one-qs2)*(
one-qs2**nstr2))
2858 dxfirst2=lenstr2*(
one-qs2)/(
one-qs2**nstr2)
2861 if(nuni2 .gt. 0)
then
2862 do ix2=nstr2+1,nstr2+nuni2
2864 xif2(ix2)=lenstr2+(dble(ix2)-0.5d0-nstr2)*dxif2(ix2)+xifmin2
2869 dxif2(ix2)=dxfirst2*qs2**(nstr2-ix2)
2870 xif2(ix2)=xifmin2+lenstr2-dxif2(ix2)*
half-dxfirst2*(
one-qs2**(nstr2-ix2))/(
one-qs2)
2873 do ix2=nstr2+nuni2+1,nxif2
2874 dxif2(ix2)=dxfirst2*qs2**(ix2-nstr2-nuni2-1)
2875 xif2(ix2)=xifmax2-lenstr2+dxif2(ix2)*
half+dxfirst2*(
one-qs2**(ix2-nstr2-nuni2-1))/(
one-qs2)
2878 call mpistop(
"unknown stretch type")
2881 if (
mype==0 .and. datatype==
'image_euv')
then
2889 write(*,
'(a,i8,a,i8)')
' Native data-resolution image grid: ',nxif1,
' x ',nxif2
2890 write(*,
'(a,f10.3,a,f10.3,a,f8.3,a,f8.3,a)') &
2891 ' Native xI1 pixel-size range: ',minval(dxif1)*length_to_km,
'--', &
2892 maxval(dxif1)*length_to_km,
' km (',minval(dxif1)/arcsec,
'--',maxval(dxif1)/arcsec,
' arcsec)'
2893 write(*,
'(a,f10.3,a,f10.3,a,f8.3,a,f8.3,a)') &
2894 ' Native xI2 pixel-size range: ',minval(dxif2)*length_to_km,
'--', &
2895 maxval(dxif2)*length_to_km,
' km (',minval(dxif2)/arcsec,
'--',maxval(dxif2)/arcsec,
' arcsec)'
2899 if (datatype==
'image_euv')
then
2908 allocate(wi(nxif1,nxif2,numwi))
2909 allocate(euv(nxif1,nxif2),dpl(nxif1,nxif2))
2911 allocate(euvthin(nxif1,nxif2),tau(nxif1,nxif2))
2920 if (has_doppler_output)
then
2925 allocate(euvs(nxif1,nxif2),dpls(nxif1,nxif2))
2934 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision, &
2939 do iigrid=1,igridstail; igrid=igrids(iigrid);
2943 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision, &
2945 call mpi_allreduce(dpls,dpl,numsi,mpi_double_precision, &
2948 if (has_doppler_output)
then
2952 deallocate(euvs,dpls)
2954 if (has_thick_output)
then
2955 if (has_doppler_output)
then
2957 has_thick_output,dpl=dpl,tau=tau,euvthin=euvthin)
2960 has_thick_output,tau=tau,euvthin=euvthin)
2962 else if (has_doppler_output)
then
2964 has_thick_output,dpl=dpl)
2974 euv,nxip1,nxip2,xip1,xip2,&
2975 dxip1,dxip2,wip,numwip,tau=tau,brightthin=euvthin)
2978 euv,nxip1,nxip2,xip1,xip2,&
2979 dxip1,dxip2,wip,numwip)
2983 euv,dpl,nxip1,nxip2,xip1,xip2,&
2984 dxip1,dxip2,wip,numwip,tau=tau,euvthin=euvthin)
2987 euv,dpl,nxip1,nxip2,xip1,xip2,&
2988 dxip1,dxip2,wip,numwip)
2990 call output_data(qunit,xip1,xip2,dxip1,dxip2,wip,nxip1,nxip2,numwip,datatype)
2991 deallocate(xip1,xip2,dxip1,dxip2,wip)
2993 call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
2996 deallocate(wi,euv,dpl,euvthin,tau)
2998 deallocate(wi,euv,dpl)
3003 if (datatype==
'image_sxr')
then
3011 allocate(wi(nxif1,nxif2,numwi))
3012 allocate(sxrs(nxif1,nxif2),sxr(nxif1,nxif2))
3015 do iigrid=1,igridstail; igrid=igrids(iigrid);
3019 call mpi_allreduce(sxrs,sxr,numsi,mpi_double_precision, &
3022 sxr=sxr*(rhessi_rsl*arcsec)**2
3030 call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
3031 deallocate(wi,sxr,sxrs)
3034 deallocate(xif1,xif2,dxif1,dxif2)
3041 integer,
intent(in) :: igrid,nXIF1,nXIF2
3042 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3043 double precision,
intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
3045 double precision,
intent(out) :: SXR(nXIF1,nXIF2)
3047 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
3048 double precision :: xb^L,xd^D
3049 double precision,
allocatable :: flux(:^D&),opacity(:^D&)
3050 double precision,
allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
3051 double precision,
allocatable :: SXRg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
3052 integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
3053 double precision :: SXRt,xc^L,xg^L,r2,area_1AU
3054 integer :: ixP^L,ixP^D
3055 integer :: direction_LOS
3065 ^d&ixomin^d=ixmlo^d\
3066 ^d&ixomax^d=ixmhi^d\
3067 ^d&iximin^d=
ixglo^d\
3068 ^d&iximax^d=
ixghi^d\
3072 allocate(flux(ixi^s))
3073 allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
3074 dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
3075 dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
3076 dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
3081 levelg=ps(igrid)%level
3085 select case(direction_los)
3096 allocate(sxrg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
3102 select case(direction_los)
3104 do ix2=ixomin2,ixomax2
3105 ixgmin1=(ix2-1)*rft+1
3107 do ix3=ixomin3,ixomax3
3108 ixgmin2=(ix3-1)*rft+1
3111 do ix1=ixomin1,ixomax1
3114 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
3118 do ix3=ixomin3,ixomax3
3119 ixgmin1=(ix3-1)*rft+1
3121 do ix1=ixomin1,ixomax1
3122 ixgmin2=(ix1-1)*rft+1
3125 do ix2=ixomin2,ixomax2
3128 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
3132 do ix1=ixomin1,ixomax1
3133 ixgmin1=(ix1-1)*rft+1
3135 do ix2=ixomin2,ixomax2
3136 ixgmin2=(ix2-1)*rft+1
3139 do ix3=ixomin3,ixomax3
3142 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
3152 select case(direction_los)
3154 ixgmin1=(ixomin2-1)*rft+1
3156 ixgmin2=(ixomin3-1)*rft+1
3159 ixgmin1=(ixomin3-1)*rft+1
3161 ixgmin2=(ixomin1-1)*rft+1
3164 ixgmin1=(ixomin1-1)*rft+1
3166 ixgmin2=(ixomin2-1)*rft+1
3170 select case(direction_los)
3172 ixpmin1=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3173 ixpmax1=
node(pig2_,igrid)*rft*block_nx2
3174 ixpmin2=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3175 ixpmax2=
node(pig3_,igrid)*rft*block_nx3
3177 ixpmin1=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3178 ixpmax1=
node(pig3_,igrid)*rft*block_nx3
3179 ixpmin2=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3180 ixpmax2=
node(pig1_,igrid)*rft*block_nx1
3182 ixpmin1=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3183 ixpmax1=
node(pig1_,igrid)*rft*block_nx1
3184 ixpmin2=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3185 ixpmax2=
node(pig2_,igrid)*rft*block_nx2
3187 xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
3188 xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
3189 dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
3190 dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
3191 sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
3192 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
3194 deallocate(flux,dxb1,dxb2,dxb3,sxrg,xg1,xg2,dxg1,dxg2)
3201 integer,
intent(in) :: igrid,nXIF1,nXIF2
3202 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3203 double precision,
intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
3205 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
3207 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
3208 double precision :: xb^L,xd^D
3209 double precision,
allocatable :: flux(:^D&),v(:^D&),rho(:^D&),opacity(:^D&)
3210 double precision,
allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
3211 double precision,
allocatable :: EUVg(:,:),Fvg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
3212 integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
3213 double precision :: EUVt,Fvt,xc^L,xg^L,r2
3214 integer :: ixP^L,ixP^D
3215 integer :: direction_LOS
3225 ^d&ixomin^d=ixmlo^d\
3226 ^d&ixomax^d=ixmhi^d\
3227 ^d&iximin^d=
ixglo^d\
3228 ^d&iximax^d=
ixghi^d\
3232 allocate(flux(ixi^s),v(ixi^s),rho(ixi^s),opacity(ixi^s))
3233 allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
3234 dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
3235 dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
3236 dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
3249 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
3250 v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
3255 levelg=ps(igrid)%level
3259 select case(direction_los)
3270 allocate(euvg(nxg1,nxg2),fvg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
3277 select case(direction_los)
3279 do ix2=ixomin2,ixomax2
3280 ixgmin1=(ix2-1)*rft+1
3282 do ix3=ixomin3,ixomax3
3283 ixgmin2=(ix3-1)*rft+1
3287 do ix1=ixomin1,ixomax1
3291 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
3292 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
3296 do ix3=ixomin3,ixomax3
3297 ixgmin1=(ix3-1)*rft+1
3299 do ix1=ixomin1,ixomax1
3300 ixgmin2=(ix1-1)*rft+1
3304 do ix2=ixomin2,ixomax2
3308 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
3309 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
3313 do ix1=ixomin1,ixomax1
3314 ixgmin1=(ix1-1)*rft+1
3316 do ix2=ixomin2,ixomax2
3317 ixgmin2=(ix2-1)*rft+1
3321 do ix3=ixomin3,ixomax3
3325 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
3326 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
3337 select case(direction_los)
3339 ixgmin1=(ixomin2-1)*rft+1
3341 ixgmin2=(ixomin3-1)*rft+1
3344 ixgmin1=(ixomin3-1)*rft+1
3346 ixgmin2=(ixomin1-1)*rft+1
3349 ixgmin1=(ixomin1-1)*rft+1
3351 ixgmin2=(ixomin2-1)*rft+1
3355 select case(direction_los)
3357 ixpmin1=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3358 ixpmax1=
node(pig2_,igrid)*rft*block_nx2
3359 ixpmin2=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3360 ixpmax2=
node(pig3_,igrid)*rft*block_nx3
3362 ixpmin1=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3363 ixpmax1=
node(pig3_,igrid)*rft*block_nx3
3364 ixpmin2=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3365 ixpmax2=
node(pig1_,igrid)*rft*block_nx1
3367 ixpmin1=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3368 ixpmax1=
node(pig1_,igrid)*rft*block_nx1
3369 ixpmin2=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3370 ixpmax2=
node(pig2_,igrid)*rft*block_nx2
3372 xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
3373 xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
3374 dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
3375 dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
3376 euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
3377 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
3378 dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
3379 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
3381 deallocate(flux,v,opacity,dxb1,dxb2,dxb3,euvg,fvg,xg1,xg2,dxg1,dxg2)
3390 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),box_min(1:3),box_max(1:3)
3391 logical,
intent(out) :: hit
3392 double precision,
intent(out) :: t_enter,t_exit
3395 double precision :: t1,t2,td
3401 if (abs(ray_dir(idir))<=smalldouble)
then
3402 if (ray_origin(idir)<box_min(idir) .or. ray_origin(idir)>box_max(idir))
then
3407 t1=(box_min(idir)-ray_origin(idir))/ray_dir(idir)
3408 t2=(box_max(idir)-ray_origin(idir))/ray_dir(idir)
3414 t_enter=max(t_enter,t1)
3415 t_exit=min(t_exit,t2)
3416 if (t_enter>=t_exit)
then
3427 integer,
intent(in) :: ixI^L, ixO^L
3428 double precision,
intent(in) :: x(ixI^S,1:ndim),dx(ixI^S,1:ndim)
3429 double precision,
allocatable,
intent(out) :: xface1(:),xface2(:),xface3(:)
3433 allocate(xface1(ixomin1:ixomax1+1),xface2(ixomin2:ixomax2+1),xface3(ixomin3:ixomax3+1))
3437 do ix1=ixomin1,ixomax1
3438 xface1(ix1)=x(ix^d,1)-half*dx(ix^d,1)
3440 xface1(ixomax1+1)=x(ixomax1,ixomin2,ixomin3,1)+half*dx(ixomax1,ixomin2,ixomin3,1)
3444 do ix2=ixomin2,ixomax2
3445 xface2(ix2)=x(ix^d,2)-half*dx(ix^d,2)
3447 xface2(ixomax2+1)=x(ixomin1,ixomax2,ixomin3,2)+half*dx(ixomin1,ixomax2,ixomin3,2)
3451 do ix3=ixomin3,ixomax3
3452 xface3(ix3)=x(ix^d,3)-half*dx(ix^d,3)
3454 xface3(ixomax3+1)=x(ixomin1,ixomin2,ixomax3,3)+half*dx(ixomin1,ixomin2,ixomax3,3)
3458 integer,
intent(in) :: imin,imax
3459 double precision,
intent(in) :: pos,faces(imin:imax+1)
3461 integer :: ilo,ihi,imid
3463 if (pos<=faces(imin))
then
3467 if (pos>=faces(imax+1))
then
3474 do while (ihi-ilo>1)
3476 if (pos>=faces(imid))
then
3482 idx=min(imax,max(imin,ilo))
3486 integer,
intent(in) :: imin,imax,idx
3487 double precision,
intent(in) :: ray_origin_axis,ray_dir_axis,faces(imin:imax+1)
3488 integer,
intent(out) :: step
3489 double precision,
intent(out) :: tMax
3491 if (ray_dir_axis>zero)
then
3493 tmax=(faces(idx+1)-ray_origin_axis)/ray_dir_axis
3494 else if (ray_dir_axis<zero)
then
3496 tmax=(faces(idx)-ray_origin_axis)/ray_dir_axis
3504 integer,
intent(in) :: imin,imax,step
3505 double precision,
intent(in) :: ray_origin_axis,ray_dir_axis,faces(imin:imax+1)
3506 integer,
intent(inout) :: idx
3507 double precision,
intent(inout) :: tMax
3508 logical,
intent(out) :: done
3512 if (idx<imin .or. idx>imax)
then
3517 tmax=(faces(idx+1)-ray_origin_axis)/ray_dir_axis
3518 else if (step<0)
then
3519 tmax=(faces(idx)-ray_origin_axis)/ray_dir_axis
3526 ray_origin,xface1,xface2,xface3,t_enter,t_exit,EUVp,Dplp)
3527 integer,
intent(in) :: ixI^L, ixO^L
3528 double precision,
intent(in) :: source(ixI^S),sourcev(ixI^S)
3529 double precision,
intent(in) :: ray_origin(1:3)
3530 double precision,
intent(in) :: xface1(ixOmin1:ixOmax1+1),xface2(ixOmin2:ixOmax2+1),&
3531 xface3(ixOmin3:ixOmax3+1)
3532 double precision,
intent(in) :: t_enter,t_exit
3533 double precision,
intent(inout) :: EUVp,Dplp
3535 integer :: ix^D,step(1:3)
3536 double precision :: pos(1:3),tMax(1:3),tNow,tNext,ds_cm,epsRay
3539 if (t_exit<=t_enter)
return
3540 epsray=max(1.d-12,1.d-10*abs(t_exit-t_enter))
3541 pos=ray_origin+(t_enter+epsray)*
vec_los
3552 tnext=min(t_exit,tmax(1),tmax(2),tmax(3))
3553 if (tnext>tnow)
then
3554 ds_cm=(tnext-tnow)*unit_length
3555 if (si_unit) ds_cm=ds_cm*1.d2
3556 euvp=euvp+source(ix^d)*ds_cm
3557 dplp=dplp+sourcev(ix^d)*ds_cm
3560 if (tnow>=t_exit-epsray)
exit
3562 if (tmax(1)<=tnow+epsray)
then
3566 if (tmax(2)<=tnow+epsray)
then
3570 if (tmax(3)<=tnow+epsray)
then
3578 double precision,
allocatable,
intent(inout) :: segments(:,:)
3579 integer,
intent(inout) :: nseg,capacity
3580 integer,
intent(in) :: pixel_id
3581 double precision,
intent(in) :: tseg,jds,kds,jvds
3583 double precision,
allocatable :: tmp(:,:)
3584 integer :: new_capacity
3586 if (capacity<=0)
then
3588 allocate(segments(5,capacity))
3589 else if (nseg>=capacity)
then
3590 new_capacity=2*capacity
3591 allocate(tmp(5,new_capacity))
3592 tmp(:,1:capacity)=segments(:,1:capacity)
3593 call move_alloc(tmp,segments)
3594 capacity=new_capacity
3598 segments(1,nseg)=dble(pixel_id)
3599 segments(2,nseg)=tseg
3600 segments(3,nseg)=jds
3601 segments(4,nseg)=kds
3602 segments(5,nseg)=jvds
3606 pixel_id,ray_origin,xface1,xface2,xface3,t_enter,t_exit,&
3607 segments,nseg,capacity)
3608 integer,
intent(in) :: ixI^L, ixO^L
3609 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S),sourcev(ixI^S)
3610 integer,
intent(in) :: pixel_id
3611 double precision,
intent(in) :: ray_origin(1:3)
3612 double precision,
intent(in) :: xface1(ixOmin1:ixOmax1+1),xface2(ixOmin2:ixOmax2+1),&
3613 xface3(ixOmin3:ixOmax3+1)
3614 double precision,
intent(in) :: t_enter,t_exit
3615 double precision,
allocatable,
intent(inout) :: segments(:,:)
3616 integer,
intent(inout) :: nseg,capacity
3618 integer :: ix^D,step(1:3)
3619 double precision :: pos(1:3),tMax(1:3),tNow,tNext,ds_cm,epsRay,tseg
3620 double precision :: jds,kds,jvds
3623 if (t_exit<=t_enter)
return
3624 epsray=max(1.d-12,1.d-10*abs(t_exit-t_enter))
3625 pos=ray_origin+(t_enter+epsray)*
vec_los
3636 tnext=min(t_exit,tmax(1),tmax(2),tmax(3))
3637 if (tnext>tnow)
then
3638 ds_cm=(tnext-tnow)*unit_length
3639 if (si_unit) ds_cm=ds_cm*1.d2
3640 jds=source(ix^d)*ds_cm
3641 kds=opacity(ix^d)*ds_cm
3642 jvds=sourcev(ix^d)*ds_cm
3643 if (jds/=zero .or. kds/=zero .or. jvds/=zero)
then
3644 tseg=half*(tnow+tnext)
3649 if (tnow>=t_exit-epsray)
exit
3651 if (tmax(1)<=tnow+epsray)
then
3655 if (tmax(2)<=tnow+epsray)
then
3659 if (tmax(3)<=tnow+epsray)
then
3667 double precision,
intent(in) :: segments(:,:)
3668 integer,
intent(inout) :: idx(:)
3669 integer,
intent(in) :: nidx
3680 double precision,
intent(in) :: segments(:,:)
3681 integer,
intent(inout) :: idx(:)
3682 integer,
intent(in) :: ilo,ihi
3689 do while (j>=ilo .and. segments(2,idx(j))>segments(2,key))
3698 double precision,
intent(in) :: segments(:,:)
3699 integer,
intent(inout) :: idx(:)
3700 integer,
intent(in) :: ilo,ihi
3703 double precision :: pivot
3705 if (ihi-ilo<=32)
then
3712 pivot=segments(2,idx((ilo+ihi)/2))
3714 do while (segments(2,idx(i))<pivot)
3717 do while (segments(2,idx(j))>pivot)
3735 integer,
intent(in) :: pixel_id
3737 owner=mod(pixel_id-1,npe)
3741 double precision,
intent(in) :: segments(:,:)
3742 integer,
intent(in) :: is,nvars
3748 if (segments(iv,is)/=segments(iv,is) .or. abs(segments(iv,is))>=1.d90)
then
3755 subroutine cart_dda_block_pixel_range(box_min,box_max,nXIF1,nXIF2,xIF1,xIF2,ixPmin1,ixPmax1,ixPmin2,ixPmax2,has_pixels)
3756 double precision,
intent(in) :: box_min(1:3),box_max(1:3)
3757 integer,
intent(in) :: nXIF1,nXIF2
3758 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3759 integer,
intent(out) :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
3760 logical,
intent(out) :: has_pixels
3763 double precision :: vec_cor(1:3),xI_cor(1:2)
3764 double precision :: xmin1,xmax1,xmin2,xmax2,dx1,dx2
3767 if (i1==1) vec_cor(1)=box_min(1)
3768 if (i1==2) vec_cor(1)=box_max(1)
3770 if (i2==1) vec_cor(2)=box_min(2)
3771 if (i2==2) vec_cor(2)=box_max(2)
3773 if (i3==1) vec_cor(3)=box_min(3)
3774 if (i3==2) vec_cor(3)=box_max(3)
3776 if (i1==1 .and. i2==1 .and. i3==1)
then
3782 xmin1=min(xmin1,xi_cor(1))
3783 xmax1=max(xmax1,xi_cor(1))
3784 xmin2=min(xmin2,xi_cor(2))
3785 xmax2=max(xmax2,xi_cor(2))
3792 dx1=abs(xif1(2)-xif1(1))
3794 dx1=max(one,abs(xmax1-xmin1))
3797 dx2=abs(xif2(2)-xif2(1))
3799 dx2=max(one,abs(xmax2-xmin2))
3802 ixpmin1=max(1,floor((xmin1-xif1(1))/dx1)+1-1)
3803 ixpmax1=min(nxif1,ceiling((xmax1-xif1(1))/dx1)+1+1)
3804 ixpmin2=max(1,floor((xmin2-xif2(1))/dx2)+1-1)
3805 ixpmax2=min(nxif2,ceiling((xmax2-xif2(1))/dx2)+1+1)
3806 has_pixels=ixpmin1<=ixpmax1 .and. ixpmin2<=ixpmax2
3812 integer,
intent(in) :: nXIF1,nXIF2
3813 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3815 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
3817 integer :: ixO^L,ixO^D,ixI^L,ix^D
3818 integer :: iigrid,igrid,ixP1,ixP2,numSI,ixPmin1,ixPmax1,ixPmin2,ixPmax2
3819 double precision :: box_min(1:3),box_max(1:3),ray_origin(1:3)
3820 double precision :: t_enter,t_exit,vlos
3821 double precision :: profile_local(2),profile_global(2)
3822 logical :: hit,has_pixels
3823 double precision,
allocatable :: source(:^D&),sourcev(:^D&),rho(:^D&),opacity(:^D&)
3824 double precision,
allocatable :: xface1(:),xface2(:),xface3(:)
3825 double precision,
allocatable :: EUVs(:,:),Dpls(:,:)
3827 allocate(euvs(nxif1,nxif2),dpls(nxif1,nxif2))
3832 do iigrid=1,igridstail; igrid=igrids(iigrid);
3833 ^d&ixomin^d=ixmlo^d\
3834 ^d&ixomax^d=ixmhi^d\
3835 ^d&iximin^d=
ixglo^d\
3836 ^d&iximax^d=
ixghi^d\
3838 box_min(1)=
rnode(rpxmin1_,igrid)
3839 box_min(2)=
rnode(rpxmin2_,igrid)
3840 box_min(3)=
rnode(rpxmin3_,igrid)
3841 box_max(1)=
rnode(rpxmax1_,igrid)
3842 box_max(2)=
rnode(rpxmax2_,igrid)
3843 box_max(3)=
rnode(rpxmax3_,igrid)
3845 xface1,xface2,xface3)
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)
4438 integer,
intent(in) :: ixI^L, ixO^L
4439 double precision,
intent(in) :: x(ixI^S,1:ndim),dx(ixI^S,1:ndim)
4440 double precision,
allocatable,
intent(out) :: rface(:),thetaface(:),phiface(:)
4441 integer :: ix1,ix2,ix3
4443 allocate(rface(ixomin1:ixomax1+1),thetaface(ixomin2:ixomax2+1),phiface(ixomin3:ixomax3+1))
4444 do ix1=ixomin1,ixomax1
4445 rface(ix1)=x(ix1,ixomin2,ixomin3,1)-half*dx(ix1,ixomin2,ixomin3,1)
4447 rface(ixomax1+1)=x(ixomax1,ixomin2,ixomin3,1)+half*dx(ixomax1,ixomin2,ixomin3,1)
4448 do ix2=ixomin2,ixomax2
4449 thetaface(ix2)=x(ixomin1,ix2,ixomin3,2)-half*dx(ixomin1,ix2,ixomin3,2)
4451 thetaface(ixomax2+1)=x(ixomin1,ixomax2,ixomin3,2)+half*dx(ixomin1,ixomax2,ixomin3,2)
4452 do ix3=ixomin3,ixomax3
4453 phiface(ix3)=x(ixomin1,ixomin2,ix3,3)-half*dx(ixomin1,ixomin2,ix3,3)
4455 phiface(ixomax3+1)=x(ixomin1,ixomin2,ixomax3,3)+half*dx(ixomin1,ixomin2,ixomax3,3)
4459 double precision,
allocatable,
intent(inout) :: tvals(:)
4460 integer,
intent(inout) :: nt,capacity
4461 double precision,
intent(in) :: t
4463 double precision,
allocatable :: tmp(:)
4465 if (t /= t .or. abs(t)>1.d90)
return
4466 if (.not.
allocated(tvals))
then
4468 allocate(tvals(capacity))
4469 else if (nt>=capacity)
then
4470 allocate(tmp(capacity))
4473 allocate(tvals(2*capacity))
4474 tvals(1:capacity)=tmp
4483 double precision,
intent(inout) :: tvals(:)
4484 integer,
intent(inout) :: nt
4485 double precision,
intent(in) :: t
4487 if (t /= t .or. abs(t)>1.d90)
return
4488 if (nt>=
size(tvals))
return
4494 double precision,
intent(inout) :: tvals(:)
4495 integer,
intent(inout) :: nt
4498 double precision :: key,epsT
4504 do while (j>=1 .and. tvals(j)>key)
4510 epst=max(1.d-12,1.d-10*max(one,abs(tvals(nt)-tvals(1))))
4513 if (abs(tvals(i)-tvals(nout))>epst)
then
4515 tvals(nout)=tvals(i)
4522 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),rface
4523 double precision,
allocatable,
intent(inout) :: tvals(:)
4524 integer,
intent(inout) :: nt,capacity
4526 double precision :: aa,bb,cc,disc,root
4529 bb=2.d0*sum(ray_origin*ray_dir)
4530 cc=sum(ray_origin**2)-rface**2
4531 disc=bb**2-4.d0*aa*cc
4532 if (disc<zero)
return
4533 root=sqrt(max(zero,disc))
4534 call sph_add_t(tvals,nt,capacity,(-bb-root)/(2.d0*aa))
4535 call sph_add_t(tvals,nt,capacity,(-bb+root)/(2.d0*aa))
4539 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),thetaface
4540 double precision,
allocatable,
intent(inout) :: tvals(:)
4541 integer,
intent(inout) :: nt,capacity
4543 double precision :: cth,aa,bb,cc,disc,root
4549 if (abs(cth)<1.d-12)
then
4550 if (abs(ray_dir(3))>1.d-14) &
4551 call sph_add_t(tvals,nt,capacity,-ray_origin(3)/ray_dir(3))
4554 aa=ray_dir(3)**2-cth**2*sum(ray_dir**2)
4555 bb=2.d0*(ray_origin(3)*ray_dir(3)-cth**2*sum(ray_origin*ray_dir))
4556 cc=ray_origin(3)**2-cth**2*sum(ray_origin**2)
4557 if (abs(aa)<1.d-14)
then
4558 if (abs(bb)>1.d-14)
call sph_add_t(tvals,nt,capacity,-cc/bb)
4561 disc=bb**2-4.d0*aa*cc
4562 if (disc<zero)
return
4563 root=sqrt(max(zero,disc))
4564 call sph_add_t(tvals,nt,capacity,(-bb-root)/(2.d0*aa))
4565 call sph_add_t(tvals,nt,capacity,(-bb+root)/(2.d0*aa))
4569 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),phiface
4570 double precision,
allocatable,
intent(inout) :: tvals(:)
4571 integer,
intent(inout) :: nt,capacity
4573 double precision :: normal(1:3),denom,numer
4575 normal(1)=-sin(phiface)
4576 normal(2)=cos(phiface)
4578 denom=sum(normal*ray_dir)
4579 if (abs(denom)<1.d-14)
return
4580 numer=sum(normal*ray_origin)
4581 call sph_add_t(tvals,nt,capacity,-numer/denom)
4585 double precision,
intent(in) :: pos(1:3)
4586 double precision,
intent(out) :: sph(1:3)
4588 sph(1)=sqrt(sum(pos**2))
4589 if (sph(1)>zero)
then
4590 sph(2)=acos(max(-one,min(one,pos(3)/sph(1))))
4594 sph(3)=atan2(pos(2),pos(1))
4598 integer,
intent(in) :: imin,imax
4599 double precision,
intent(in) ::
value,faces(imin:imax+1)
4601 integer :: ilo,ihi,imid
4604 if (
value<faces(imin)-1.d-12 .or.
value>faces(imax+1)+1.d-12)
return
4605 if (
value<=faces(imin))
then
4609 if (
value>=faces(imax+1))
then
4616 do while (ihi-ilo>1)
4618 if (
value>=faces(imid))
then
4624 idx=min(imax,max(imin,ilo))
4628 integer,
intent(in) :: imin,imax
4629 double precision,
intent(in) ::
value,faces(imin:imax+1)
4631 integer :: ilo,ihi,imid
4634 if (
value>faces(imin)+1.d-12 .or.
value<faces(imax+1)-1.d-12)
return
4635 if (
value>=faces(imin))
then
4639 if (
value<=faces(imax+1))
then
4646 do while (ihi-ilo>1)
4648 if (
value<=faces(imid))
then
4654 idx=min(imax,max(imin,ilo))
4658 double precision,
intent(in) :: pos(1:3)
4659 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4660 phiface(ixOmin3:ixOmax3+1)
4661 integer,
intent(in) :: ixO^L
4662 integer,
intent(out) :: ix1,ix2,ix3
4663 logical,
intent(out) :: inside
4665 double precision :: sph(1:3),phi
4669 if (phi<phiface(ixomin3)-1.d-12) phi=phi+2.d0*dpi
4670 if (phi>phiface(ixomax3+1)+1.d-12) phi=phi-2.d0*dpi
4674 inside=ix1>0 .and. ix2>0 .and. ix3>0
4679 double precision,
intent(in) :: pos(1:3),ximg1,ximg2
4681 double precision :: dotp,rc,rthick,rloc
4684 rc=sqrt(ximg1**2+ximg2**2)
4685 rloc=sqrt(sum(pos**2))
4688 if (dotp>=
zero)
then
4689 if (rc<=rthick) visible=.false.
4691 if (rloc<=rthick) visible=.false.
4696 ixPmin1,ixPmax1,ixPmin2,ixPmax2,has_pixels)
4697 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4698 phiface(ixOmin3:ixOmax3+1)
4699 integer,
intent(in) :: ixO^L,nXI1,nXI2
4700 double precision,
intent(in) :: xI1(nXI1),xI2(nXI2),dxI
4701 integer,
intent(out) :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
4702 logical,
intent(out) :: has_pixels
4704 integer,
parameter :: nsample=5
4706 double precision :: sph(1:3),xcent(1:2)
4707 double precision :: xmin1,xmax1,xmin2,xmax2
4708 double precision :: wr,wt,wp,pad
4716 wr=dble(ir)/dble(nsample-1)
4717 sph(1)=(one-wr)*rface(ixomin1)+wr*rface(ixomax1+1)
4719 wt=dble(it)/dble(nsample-1)
4720 sph(2)=(one-wt)*thetaface(ixomin2)+wt*thetaface(ixomax2+1)
4722 wp=dble(ip)/dble(nsample-1)
4723 if (ir/=0 .and. ir/=nsample-1 .and. it/=0 .and. it/=nsample-1 .and. &
4724 ip/=0 .and. ip/=nsample-1) cycle
4725 sph(3)=(one-wp)*phiface(ixomin3)+wp*phiface(ixomax3+1)
4727 xmin1=min(xmin1,xcent(1))
4728 xmax1=max(xmax1,xcent(1))
4729 xmin2=min(xmin2,xcent(2))
4730 xmax2=max(xmax2,xcent(2))
4739 ixpmin1=max(1,floor((xmin1-(xi1(1)-half*dxi))/dxi)+1)
4740 ixpmax1=min(nxi1,ceiling((xmax1-(xi1(1)-half*dxi))/dxi))
4741 ixpmin2=max(1,floor((xmin2-(xi2(1)-half*dxi))/dxi)+1)
4742 ixpmax2=min(nxi2,ceiling((xmax2-(xi2(1)-half*dxi))/dxi))
4743 has_pixels=ixpmin1<=ixpmax1 .and. ixpmin2<=ixpmax2
4747 rface,thetaface,phiface,EUVp)
4748 integer,
intent(in) :: ixI^L,ixO^L
4749 double precision,
intent(in) :: source(ixI^S),ray_origin(1:3),ximg1,ximg2
4750 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4751 phiface(ixOmin3:ixOmax3+1)
4752 double precision,
intent(inout) :: EUVp
4754 integer :: nt,capacity,i,ix^D
4755 double precision,
allocatable :: tvals(:)
4756 double precision :: posMid(1:3),ds_cm,tMid,t0,t1
4761 do ix1=ixomin1,ixomax1+1
4764 do ix2=ixomin2,ixomax2+1
4767 do ix3=ixomin3,ixomax3+1
4771 if (
allocated(tvals))
deallocate(tvals)
4781 posmid=ray_origin+tmid*
vec_los
4783 call sph_locate_cell(posmid,rface,thetaface,phiface,ixo^l,ix1,ix2,ix3,inside)
4784 if (.not. inside) cycle
4785 ds_cm=(t1-t0)*unit_length
4786 if (si_unit) ds_cm=ds_cm*1.d2
4787 euvp=euvp+source(ix^d)*ds_cm
4793 pixel_id,ray_origin,ximg1,ximg2,&
4794 rface,thetaface,phiface,rface2,&
4795 theta_cos,phi_sin,phi_cos,&
4796 segments,nseg,capacity)
4799 integer,
intent(in) :: ixI^L,ixO^L,pixel_id
4800 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S)
4801 double precision,
intent(in) :: ray_origin(1:3),ximg1,ximg2
4802 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4803 phiface(ixOmin3:ixOmax3+1)
4804 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4805 phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
4806 double precision,
allocatable,
intent(inout) :: segments(:,:)
4807 integer,
intent(inout) :: nseg,capacity
4809 integer :: nt,i,ix^D
4810 double precision :: tvals(2*(ixOmax1-ixOmin1+2)+2*(ixOmax2-ixOmin2+2)+&
4811 (ixOmax3-ixOmin3+2))
4812 double precision :: posMid(1:3),ds_cm,tMid,t0,t1,jds,kds
4813 double precision :: dir2,origin2,odotd,aa,bb,cc,disc,root,cth2
4814 double precision :: denom,numer,r2,mu,phi,dotp,rthick2,rc2
4819 origin2=sum(ray_origin**2)
4821 rthick2=(r_opt_thick*
const_rsun/unit_length)**2
4822 rc2=ximg1**2+ximg2**2
4824 do ix1=ixomin1,ixomax1+1
4827 cc=origin2-rface2(ix1)
4828 disc=bb**2-4.d0*aa*cc
4829 if (disc>=
zero)
then
4830 root=sqrt(max(
zero,disc))
4835 do ix2=ixomin2,ixomax2+1
4836 if (abs(theta_cos(ix2))<1.d-12)
then
4841 cth2=theta_cos(ix2)**2
4843 bb=2.d0*(ray_origin(3)*
vec_los(3)-cth2*odotd)
4844 cc=ray_origin(3)**2-cth2*origin2
4845 if (abs(aa)<1.d-14)
then
4848 disc=bb**2-4.d0*aa*cc
4849 if (disc>=
zero)
then
4850 root=sqrt(max(
zero,disc))
4856 do ix3=ixomin3,ixomax3+1
4858 if (abs(denom)>=1.d-14)
then
4859 numer=-phi_sin(ix3)*ray_origin(1)+phi_cos(ix3)*ray_origin(2)
4871 posmid=ray_origin+tmid*
vec_los
4874 if (dotp>=
zero)
then
4875 if (rc2<=rthick2) cycle
4877 if (r2<=rthick2) cycle
4882 mu=posmid(3)/sqrt(r2)
4887 phi=atan2(posmid(2),posmid(1))
4888 if (phi<phiface(ixomin3)-1.d-12) phi=phi+2.d0*
dpi
4889 if (phi>phiface(ixomax3+1)+1.d-12) phi=phi-2.d0*
dpi
4891 inside=ix1>0 .and. ix2>0 .and. ix3>0
4892 if (.not. inside) cycle
4893 ds_cm=(t1-t0)*unit_length
4894 if (si_unit) ds_cm=ds_cm*1.d2
4895 jds=max(
zero,source(ix^d))*ds_cm
4896 kds=max(
zero,opacity(ix^d))*ds_cm
4902 double precision,
intent(in) :: pos(1:3)
4903 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4904 phiface(ixOmin3:ixOmax3+1)
4905 integer,
intent(in) :: ixO^L
4906 integer,
intent(out) :: ix1,ix2,ix3
4907 logical,
intent(out) :: inside
4909 double precision :: r2,mu,phi
4919 phi=atan2(pos(2),pos(1))
4920 if (phi<phiface(ixomin3)-1.d-12) phi=phi+2.d0*dpi
4921 if (phi>phiface(ixomax3+1)+1.d-12) phi=phi-2.d0*dpi
4923 inside=ix1>0 .and. ix2>0 .and. ix3>0
4927 double precision,
intent(in) :: t,tNow,tExit,epsRay
4928 double precision,
intent(inout) :: tNext
4929 logical,
intent(inout) :: found
4931 if (t>tnow+epsray .and. t<=texit+epsray .and. t<tnext)
then
4938 double precision,
intent(in) :: t,theta_face_cos,ray_origin(1:3),tNow,tExit,epsRay
4939 double precision,
intent(inout) :: tNext
4940 logical,
intent(inout) :: found
4942 double precision :: pos(1:3),r2
4944 if (t<=tnow+epsray .or. t>texit+epsray .or. t>=tnext)
return
4947 if (r2<=zero)
return
4948 if (theta_face_cos>1.d-12 .and. pos(3)<-1.d-10)
return
4949 if (theta_face_cos<-1.d-12 .and. pos(3)>1.d-10)
return
4950 if (abs(pos(3)**2-theta_face_cos**2*r2)>1.d-6*max(one,r2))
return
4956 double precision,
intent(in) :: t,phi_face_sin,phi_face_cos,ray_origin(1:3),tNow,tExit,epsRay
4957 double precision,
intent(inout) :: tNext
4958 logical,
intent(inout) :: found
4960 double precision :: pos(1:3),radialDot
4962 if (t<=tnow+epsray .or. t>texit+epsray .or. t>=tnext)
return
4964 radialdot=phi_face_cos*pos(1)+phi_face_sin*pos(2)
4965 if (radialdot<-1.d-10)
return
4971 ix1,ix2,ix3,tNow,tExit,epsRay,tNext,found)
4972 double precision,
intent(in) :: ray_origin(1:3)
4973 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4974 phiface(ixOmin3:ixOmax3+1)
4975 double precision,
intent(in) :: phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
4976 integer,
intent(in) :: ixO^L,ix1,ix2,ix3
4977 double precision,
intent(in) :: tNow,tExit,epsRay
4978 double precision,
intent(out) :: tNext
4979 logical,
intent(out) :: found
4982 double precision :: dir2,origin2,odotd,aa,bb,cc,disc,root,cth2,denom,numer
4987 origin2=sum(ray_origin**2)
4993 cc=origin2-rface2(iface)
4994 disc=bb**2-4.d0*aa*cc
4995 if (disc>=zero)
then
4996 root=sqrt(max(zero,disc))
5003 if (abs(theta_cos(iface))<1.d-12)
then
5005 -ray_origin(3)/
vec_los(3),theta_cos(iface),ray_origin,tnow,texit,epsray,&
5009 cth2=theta_cos(iface)**2
5011 bb=2.d0*(ray_origin(3)*
vec_los(3)-cth2*odotd)
5012 cc=ray_origin(3)**2-cth2*origin2
5013 if (abs(aa)<1.d-14)
then
5015 ray_origin,tnow,texit,epsray,tnext,found)
5017 disc=bb**2-4.d0*aa*cc
5018 if (disc>=zero)
then
5019 root=sqrt(max(zero,disc))
5021 ray_origin,tnow,texit,epsray,tnext,found)
5023 ray_origin,tnow,texit,epsray,tnext,found)
5030 if (abs(denom)>=1.d-14)
then
5031 numer=-phi_sin(iface)*ray_origin(1)+phi_cos(iface)*ray_origin(2)
5033 tnow,texit,epsray,tnext,found)
5039 ray_origin,ximg1,ximg2,rface2,theta_cos,phiface,&
5041 t_enter,t_exit,segments,nseg,capacity,ok)
5044 integer,
intent(in) :: ixI^L,ixO^L,pixel_id
5045 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S)
5046 double precision,
intent(in) :: ray_origin(1:3),ximg1,ximg2
5047 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
5048 phiface(ixOmin3:ixOmax3+1)
5049 double precision,
intent(in) :: phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
5050 double precision,
intent(in) :: t_enter,t_exit
5051 double precision,
allocatable,
intent(inout) :: segments(:,:)
5052 integer,
intent(inout) :: nseg,capacity
5053 logical,
intent(out) :: ok
5055 integer :: ix^D,nstep,maxSteps
5056 double precision :: tNow,tNext,tEnd,tMid,epsRay,ds_cm,jds,kds
5057 double precision :: pos(1:3),r2,dotp,rthick2,rc2
5058 logical :: inside,found
5061 if (t_exit<=t_enter)
return
5062 epsray=max(1.d-12,1.d-10*max(
one,abs(t_exit-t_enter)))
5063 pos=ray_origin+(t_enter+epsray)*
vec_los
5065 if (.not. inside)
then
5070 rthick2=(r_opt_thick*
const_rsun/unit_length)**2
5071 rc2=ximg1**2+ximg2**2
5074 maxsteps=8*((ixomax1-ixomin1+1)+(ixomax2-ixomin2+1)+(ixomax3-ixomin3+1)+3)
5076 do while (tnow<t_exit-epsray)
5078 ix1,ix2,ix3,tnow,t_exit,epsray,tnext,found)
5079 if (.not. found)
then
5083 tend=min(tnext,t_exit)
5085 tmid=
half*(tnow+tend)
5089 if (.not. ((dotp>=
zero .and. rc2<=rthick2) .or. &
5090 (dotp<
zero .and. r2<=rthick2)))
then
5091 ds_cm=(tend-tnow)*unit_length
5092 if (si_unit) ds_cm=ds_cm*1.d2
5093 jds=max(
zero,source(ix^d))*ds_cm
5094 kds=max(
zero,opacity(ix^d))*ds_cm
5099 if (tnow>=t_exit-epsray)
exit
5100 pos=ray_origin+(tnow+epsray)*
vec_los
5102 if (.not. inside)
then
5107 if (nstep>maxsteps)
then
5115 pixel_id,ray_origin,ximg1,ximg2,&
5116 rface,thetaface,phiface,rface2,&
5117 theta_cos,phi_sin,phi_cos,&
5118 segments,nseg,capacity,fallback)
5119 integer,
intent(in) :: ixI^L,ixO^L,pixel_id
5120 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S)
5121 double precision,
intent(in) :: ray_origin(1:3),ximg1,ximg2
5122 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
5123 phiface(ixOmin3:ixOmax3+1)
5124 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
5125 phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
5126 double precision,
allocatable,
intent(inout) :: segments(:,:)
5127 integer,
intent(inout) :: nseg,capacity
5128 logical,
intent(out) :: fallback
5130 integer :: nt,i,nsegStart,ix^D
5131 double precision :: tvals(12),posMid(1:3),t0,t1,tMid
5132 double precision :: dir2,origin2,odotd,aa,bb,cc,disc,root,cth2,denom,numer
5133 logical :: inside,ok
5139 origin2=sum(ray_origin**2)
5142 do ix1=ixomin1,ixomax1+1,ixomax1-ixomin1+1
5145 cc=origin2-rface2(ix1)
5146 disc=bb**2-4.d0*aa*cc
5147 if (disc>=zero)
then
5148 root=sqrt(max(zero,disc))
5153 do ix2=ixomin2,ixomax2+1,ixomax2-ixomin2+1
5154 if (abs(theta_cos(ix2))<1.d-12)
then
5159 cth2=theta_cos(ix2)**2
5161 bb=2.d0*(ray_origin(3)*
vec_los(3)-cth2*odotd)
5162 cc=ray_origin(3)**2-cth2*origin2
5163 if (abs(aa)<1.d-14)
then
5166 disc=bb**2-4.d0*aa*cc
5167 if (disc>=zero)
then
5168 root=sqrt(max(zero,disc))
5174 do ix3=ixomin3,ixomax3+1,ixomax3-ixomin3+1
5176 if (abs(denom)>=1.d-14)
then
5177 numer=-phi_sin(ix3)*ray_origin(1)+phi_cos(ix3)*ray_origin(2)
5190 posmid=ray_origin+tmid*
vec_los
5192 if (.not. inside) cycle
5194 ray_origin,ximg1,ximg2,rface2,theta_cos,phiface,phi_sin,phi_cos,&
5195 t0,t1,segments,nseg,capacity,ok)
5200 pixel_id,ray_origin,ximg1,ximg2,rface,thetaface,phiface,rface2,&
5201 theta_cos,phi_sin,phi_cos,segments,nseg,capacity)
5211 integer,
intent(in) :: numXI1,numXI2
5212 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2),dxI
5214 double precision,
intent(inout) :: EM(numXI1,numXI2)
5216 integer :: ixO^L,ixI^L,ix^D
5217 integer :: iigrid,igrid,ixP1,ixP2,ixPmin1,ixPmax1,ixPmin2,ixPmax2
5218 integer :: iseg,nseg,capacity,sphDdaFallbackLocal,sphDdaFallbackGlobal
5219 double precision :: ray_origin(1:3),profile_local(3),profile_global(3)
5220 double precision,
allocatable :: source(:^D&)
5221 double precision,
allocatable :: rface(:),thetaface(:),phiface(:)
5222 double precision,
allocatable :: rface2(:),theta_cos(:),phi_sin(:),phi_cos(:)
5223 double precision,
allocatable :: segments(:,:)
5224 logical :: has_pixels,ddaFallback
5227 sphddafallbacklocal=0
5228 do iigrid=1,igridstail; igrid=igrids(iigrid);
5229 ^d&ixomin^d=ixmlo^d\
5230 ^d&ixomax^d=ixmhi^d\
5231 ^d&iximin^d=
ixglo^d\
5232 ^d&iximax^d=
ixghi^d\
5234 allocate(source(ixi^s))
5238 rface,thetaface,phiface)
5240 allocate(rface2(ixomin1:ixomax1+1),theta_cos(ixomin2:ixomax2+1),&
5241 phi_sin(ixomin3:ixomax3+1),phi_cos(ixomin3:ixomax3+1))
5243 theta_cos=cos(thetaface)
5244 phi_sin=sin(phiface)
5245 phi_cos=cos(phiface)
5248 ixpmin1,ixpmax1,ixpmin2,ixpmax2,has_pixels)
5249 if (has_pixels)
then
5251 do ixp1=ixpmin1,ixpmax1
5252 do ixp2=ixpmin2,ixpmax2
5254 profile_local(1)=profile_local(1)+one
5258 1,ray_origin,xi1(ixp1),xi2(ixp2),rface,thetaface,phiface,&
5259 rface2,theta_cos,phi_sin,phi_cos,segments,nseg,capacity,ddafallback)
5260 if (ddafallback) sphddafallbacklocal=sphddafallbacklocal+1
5262 em(ixp1,ixp2)=em(ixp1,ixp2)+segments(3,iseg)
5266 rface,thetaface,phiface,em(ixp1,ixp2))
5270 profile_local(2)=profile_local(2)+dble((ixpmax1-ixpmin1+1)*(ixpmax2-ixpmin2+1))
5272 if (
allocated(segments))
deallocate(segments)
5273 profile_local(3)=profile_local(3)+one
5274 if (
allocated(rface2))
deallocate(rface2)
5275 if (
allocated(theta_cos))
deallocate(theta_cos)
5276 if (
allocated(phi_sin))
deallocate(phi_sin)
5277 if (
allocated(phi_cos))
deallocate(phi_cos)
5278 deallocate(source,rface,thetaface,phiface)
5280 call mpi_allreduce(profile_local,profile_global,3,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5281 call mpi_allreduce(sphddafallbacklocal,sphddafallbackglobal,1,mpi_integer,mpi_sum,
icomm,
ierrmpi)
5284 write(*,
'(a,3(es12.5,1x))')
' sph_dda thin profile rays pixels blocks: ',profile_global
5285 write(*,
'(a,i0)')
' sph_dda thin fallback rays: ',sphddafallbackglobal
5287 write(*,
'(a,3(es12.5,1x))')
' sph_intersection thin profile rays pixels blocks: ',profile_global
5295 integer,
intent(in) :: numXI1,numXI2
5296 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2),dxI
5298 double precision,
intent(out) :: EUV(numXI1,numXI2),Tau(numXI1,numXI2),EUVthin(numXI1,numXI2)
5300 integer,
parameter :: nSegVars=5
5301 integer :: ixO^L,ixI^L,ix^D
5302 integer :: iigrid,igrid,ixP1,ixP2,ipix,ipixStart,ipixEnd,nPixBatch,pixel_id
5303 integer :: nseg,capacity,totalCount,totalSeg,ipe,is,iseg,nidx,owner,isegDest,nsegBefore
5304 integer :: ixGlobal,iyGlobal,ixPmin1,ixPmax1,ixPmin2,ixPmax2,iFirst,iLast,iLocal
5305 integer :: nPixBatchTarget
5306 integer :: maxSegBatchTarget,maxSegCommTarget,maxNsegBatch,nPixTotal
5307 integer :: maxOwnerSegCount,maxOwnerSegCountLocal,segOffset,recvFill,totalRoundCount,totalRoundSeg
5308 integer :: sphDdaFallbackLocal,sphDdaFallbackGlobal
5309 integer,
allocatable :: sendCounts(:),recvCounts(:),sendDispls(:),recvDispls(:)
5310 integer,
allocatable :: roundSendCounts(:),roundRecvCounts(:)
5311 integer,
allocatable :: roundSendDispls(:),roundRecvDispls(:)
5312 integer,
allocatable :: ownerSegCounts(:),ownerOffsets(:),idx(:)
5313 integer,
allocatable :: bucketCounts(:),bucketOffsets(:),bucketFill(:)
5314 double precision :: ray_origin(1:3),atten
5315 double precision :: profile_local(5),profile_global(5),profile_batch(5)
5316 double precision :: phys_max_local(2),phys_max_global(2)
5317 double precision :: phys_sum_local(2),phys_sum_global(2),phys_sum_batch(2)
5318 double precision,
allocatable :: segments(:,:),segments_send(:,:),segments_recv(:,:)
5319 double precision,
allocatable :: segments_recv_round(:,:)
5320 double precision,
allocatable :: image_reduce(:,:)
5321 logical :: has_pixels,batchAccepted,batchReduced,ddaFallback
5330 sphddafallbacklocal=0
5331 allocate(sendcounts(0:
npe-1),recvcounts(0:
npe-1),senddispls(0:
npe-1),recvdispls(0:
npe-1))
5332 allocate(roundsendcounts(0:
npe-1),roundrecvcounts(0:
npe-1))
5333 allocate(roundsenddispls(0:
npe-1),roundrecvdispls(0:
npe-1))
5334 allocate(ownersegcounts(0:
npe-1),owneroffsets(0:
npe-1))
5335 allocate(cache(igridstail))
5337 allocate(bucketcounts(npixbatchtarget),bucketoffsets(npixbatchtarget+1),&
5338 bucketfill(npixbatchtarget))
5340 do iigrid=1,igridstail; igrid=igrids(iigrid);
5341 ^d&ixomin^d=ixmlo^d\
5342 ^d&ixomax^d=ixmhi^d\
5343 ^d&iximin^d=
ixglo^d\
5344 ^d&iximax^d=
ixghi^d\
5346 cache(iigrid)%igrid=igrid
5347 allocate(cache(iigrid)%source(ixi^s),cache(iigrid)%opacity(ixi^s))
5348 cache(iigrid)%source=zero
5349 cache(iigrid)%opacity=zero
5350 call get_euv(
wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,cache(iigrid)%source)
5353 phys_max_local(1)=max(phys_max_local(1),maxval(cache(iigrid)%source(ixo^s)))
5354 phys_max_local(2)=max(phys_max_local(2),maxval(cache(iigrid)%opacity(ixo^s)))
5356 cache(iigrid)%rface,cache(iigrid)%thetaface,&
5357 cache(iigrid)%phiface)
5358 allocate(cache(iigrid)%rface2(ixomin1:ixomax1+1),&
5359 cache(iigrid)%theta_cos(ixomin2:ixomax2+1),&
5360 cache(iigrid)%phi_sin(ixomin3:ixomax3+1),&
5361 cache(iigrid)%phi_cos(ixomin3:ixomax3+1))
5362 cache(iigrid)%rface2=cache(iigrid)%rface**2
5363 cache(iigrid)%theta_cos=cos(cache(iigrid)%thetaface)
5364 cache(iigrid)%phi_sin=sin(cache(iigrid)%phiface)
5365 cache(iigrid)%phi_cos=cos(cache(iigrid)%phiface)
5367 cache(iigrid)%phiface,ixo^l,numxi1,numxi2,xi1,xi2,dxi,&
5368 cache(iigrid)%ixPmin1,cache(iigrid)%ixPmax1,&
5369 cache(iigrid)%ixPmin2,cache(iigrid)%ixPmax2,cache(iigrid)%has_pixels)
5372 npixtotal=numxi1*numxi2
5374 do while (ipixstart<=npixtotal)
5375 ipixend=min(numxi1*numxi2,ipixstart+npixbatchtarget-1)
5376 npixbatch=ipixend-ipixstart+1
5377 batchaccepted=.false.
5378 batchreduced=.false.
5380 do while (.not. batchaccepted)
5386 do iigrid=1,igridstail; igrid=igrids(iigrid);
5387 ^d&ixomin^d=ixmlo^d\
5388 ^d&ixomax^d=ixmhi^d\
5389 ^d&iximin^d=
ixglo^d\
5390 ^d&iximax^d=
ixghi^d\
5392 ixpmin1=cache(iigrid)%ixPmin1
5393 ixpmax1=cache(iigrid)%ixPmax1
5394 ixpmin2=cache(iigrid)%ixPmin2
5395 ixpmax2=cache(iigrid)%ixPmax2
5396 has_pixels=cache(iigrid)%has_pixels
5397 if (.not. has_pixels) cycle
5399 do ixp2=ixpmin2,ixpmax2
5400 ifirst=max(ipixstart,(ixp2-1)*numxi1+ixpmin1)
5401 ilast=min(ipixend,(ixp2-1)*numxi1+ixpmax1)
5402 if (ifirst>ilast) cycle
5403 do ipix=ifirst,ilast
5404 ixp1=1+mod(ipix-1,numxi1)
5407 profile_batch(1)=profile_batch(1)+one
5411 cache(iigrid)%opacity,pixel_id,ray_origin,xi1(ixp1),xi2(ixp2),&
5412 cache(iigrid)%rface,cache(iigrid)%thetaface,cache(iigrid)%phiface,&
5413 cache(iigrid)%rface2,cache(iigrid)%theta_cos,&
5414 cache(iigrid)%phi_sin,cache(iigrid)%phi_cos,&
5415 segments,nseg,capacity,ddafallback)
5416 if (ddafallback) sphddafallbacklocal=sphddafallbacklocal+1
5419 cache(iigrid)%opacity,pixel_id,ray_origin,xi1(ixp1),xi2(ixp2),&
5420 cache(iigrid)%rface,cache(iigrid)%thetaface,cache(iigrid)%phiface,&
5421 cache(iigrid)%rface2,cache(iigrid)%theta_cos,&
5422 cache(iigrid)%phi_sin,cache(iigrid)%phi_cos,&
5423 segments,nseg,capacity)
5425 if (nseg>nsegbefore) profile_batch(2)=profile_batch(2)+one
5426 profile_batch(3)=profile_batch(3)+dble(nseg-nsegbefore)
5427 do iseg=nsegbefore+1,nseg
5428 phys_sum_batch(1)=phys_sum_batch(1)+segments(3,iseg)
5429 phys_sum_batch(2)=phys_sum_batch(2)+segments(4,iseg)
5435 call mpi_allreduce(nseg,maxnsegbatch,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
5436 if (maxnsegbatch>maxsegbatchtarget .and. npixbatch>1)
then
5437 npixbatch=max(1,npixbatch/2)
5438 ipixend=ipixstart+npixbatch-1
5439 if (
allocated(segments))
deallocate(segments)
5442 batchaccepted=.true.
5446 profile_local=profile_local+profile_batch
5447 phys_sum_local=phys_sum_local+phys_sum_batch
5450 write(*,
'(a,3(i0,1x))')
' sph_dda thick adaptive batch: ',&
5451 ipixstart,ipixend,maxnsegbatch
5453 write(*,
'(a,3(i0,1x))')
' sph_intersection thick adaptive batch: ',&
5454 ipixstart,ipixend,maxnsegbatch
5458 if (.not.
allocated(segments))
then
5460 allocate(segments(nsegvars,capacity))
5465 ownersegcounts(owner)=ownersegcounts(owner)+1
5467 sendcounts=nsegvars*ownersegcounts
5470 senddispls(ipe)=senddispls(ipe-1)+sendcounts(ipe-1)
5473 allocate(segments_send(nsegvars,max(1,nseg)))
5477 isegdest=senddispls(owner)/nsegvars+owneroffsets(owner)+1
5478 segments_send(:,isegdest)=segments(:,is)
5479 owneroffsets(owner)=owneroffsets(owner)+1
5482 call mpi_alltoall(sendcounts,1,mpi_integer,recvcounts,1,mpi_integer,
icomm,
ierrmpi)
5485 recvdispls(ipe)=recvdispls(ipe-1)+recvcounts(ipe-1)
5487 totalcount=sum(recvcounts)
5488 totalseg=totalcount/nsegvars
5489 profile_local(4)=profile_local(4)+dble(totalcount)
5490 allocate(segments_recv(nsegvars,max(1,totalseg)))
5493 maxownersegcountlocal=maxval(ownersegcounts)
5494 call mpi_allreduce(maxownersegcountlocal,maxownersegcount,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
5495 do segoffset=0,maxownersegcount-1,maxsegcommtarget
5497 roundsenddispls=senddispls
5499 if (ownersegcounts(ipe)>segoffset)
then
5500 roundsendcounts(ipe)=nsegvars*min(maxsegcommtarget,ownersegcounts(ipe)-segoffset)
5501 roundsenddispls(ipe)=senddispls(ipe)+nsegvars*segoffset
5505 call mpi_alltoall(roundsendcounts,1,mpi_integer,roundrecvcounts,1,mpi_integer,
icomm,
ierrmpi)
5506 roundrecvdispls(0)=0
5508 roundrecvdispls(ipe)=roundrecvdispls(ipe-1)+roundrecvcounts(ipe-1)
5510 totalroundcount=sum(roundrecvcounts)
5511 totalroundseg=totalroundcount/nsegvars
5512 allocate(segments_recv_round(nsegvars,max(1,totalroundseg)))
5514 call mpi_alltoallv(segments_send,roundsendcounts,roundsenddispls,mpi_double_precision,&
5515 segments_recv_round,roundrecvcounts,roundrecvdispls,&
5518 if (totalroundseg>0)
then
5519 segments_recv(:,recvfill+1:recvfill+totalroundseg)=segments_recv_round(:,1:totalroundseg)
5520 recvfill=recvfill+totalroundseg
5522 deallocate(segments_recv_round)
5525 if (recvfill/=totalseg)
call mpistop(
"ray-segment receive mismatch")
5527 if (totalseg>0)
then
5528 allocate(idx(totalseg))
5529 bucketcounts(1:npixbatch)=0
5532 ipix=nint(segments_recv(1,is))
5534 ilocal=ipix-ipixstart+1
5535 bucketcounts(ilocal)=bucketcounts(ilocal)+1
5541 do ilocal=1,npixbatch
5542 bucketoffsets(ilocal+1)=bucketoffsets(ilocal)+bucketcounts(ilocal)
5544 bucketfill(1:npixbatch)=bucketoffsets(1:npixbatch)
5547 ipix=nint(segments_recv(1,is))
5549 ilocal=ipix-ipixstart+1
5550 idx(bucketfill(ilocal))=is
5551 bucketfill(ilocal)=bucketfill(ilocal)+1
5556 do ipix=ipixstart,ipixend
5558 ilocal=ipix-ipixstart+1
5559 nidx=bucketcounts(ilocal)
5561 profile_local(5)=profile_local(5)+dble(nidx)*dble(nidx)
5563 idx(bucketoffsets(ilocal):bucketoffsets(ilocal+1)-1),nidx)
5564 ixglobal=1+mod(ipix-1,numxi1)
5565 iyglobal=1+(ipix-1)/numxi1
5566 do iseg=bucketoffsets(ilocal),bucketoffsets(ilocal+1)-1
5568 euvthin(ixglobal,iyglobal)=euvthin(ixglobal,iyglobal)+segments_recv(3,is)
5570 euv(ixglobal,iyglobal)=euv(ixglobal,iyglobal)+atten*segments_recv(3,is)
5571 tau(ixglobal,iyglobal)=tau(ixglobal,iyglobal)+max(zero,segments_recv(4,is))
5578 deallocate(segments_send,segments_recv)
5579 if (
allocated(segments))
deallocate(segments)
5583 do iigrid=1,igridstail
5584 if (
allocated(cache(iigrid)%source))
deallocate(cache(iigrid)%source)
5585 if (
allocated(cache(iigrid)%opacity))
deallocate(cache(iigrid)%opacity)
5586 if (
allocated(cache(iigrid)%rface))
deallocate(cache(iigrid)%rface)
5587 if (
allocated(cache(iigrid)%thetaface))
deallocate(cache(iigrid)%thetaface)
5588 if (
allocated(cache(iigrid)%phiface))
deallocate(cache(iigrid)%phiface)
5589 if (
allocated(cache(iigrid)%rface2))
deallocate(cache(iigrid)%rface2)
5590 if (
allocated(cache(iigrid)%theta_cos))
deallocate(cache(iigrid)%theta_cos)
5591 if (
allocated(cache(iigrid)%phi_sin))
deallocate(cache(iigrid)%phi_sin)
5592 if (
allocated(cache(iigrid)%phi_cos))
deallocate(cache(iigrid)%phi_cos)
5595 deallocate(sendcounts,recvcounts,senddispls,recvdispls,roundsendcounts,roundrecvcounts,&
5596 roundsenddispls,roundrecvdispls,ownersegcounts,owneroffsets,bucketcounts,&
5597 bucketoffsets,bucketfill)
5598 allocate(image_reduce(numxi1,numxi2))
5599 call mpi_allreduce(euv,image_reduce,numxi1*numxi2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5601 call mpi_allreduce(tau,image_reduce,numxi1*numxi2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5603 call mpi_allreduce(euvthin,image_reduce,numxi1*numxi2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5604 euvthin=image_reduce
5605 deallocate(image_reduce)
5606 call mpi_allreduce(profile_local,profile_global,5,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5607 call mpi_allreduce(phys_max_local,phys_max_global,2,mpi_double_precision,mpi_max,
icomm,
ierrmpi)
5608 call mpi_allreduce(phys_sum_local,phys_sum_global,2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5609 call mpi_allreduce(sphddafallbacklocal,sphddafallbackglobal,1,mpi_integer,mpi_sum,
icomm,
ierrmpi)
5612 write(*,
'(a,5(es12.5,1x))')
' sph_dda thick profile: ',profile_global
5613 write(*,
'(a,i0)')
' sph_dda thick fallback rays: ',sphddafallbackglobal
5614 write(*,
'(a,4(es12.5,1x))')
' sph_dda thick physics maxj maxk sumjds sumkds: ',&
5615 phys_max_global(1),phys_max_global(2),phys_sum_global(1),phys_sum_global(2)
5617 write(*,
'(a,5(es12.5,1x))')
' sph_intersection thick profile: ',profile_global
5618 write(*,
'(a,4(es12.5,1x))')
' sph_intersection thick physics maxj maxk sumjds sumkds: ',&
5619 phys_max_global(1),phys_max_global(2),phys_sum_global(1),phys_sum_global(2)
5629 double precision,
intent(out) :: xImin1,xImax1,xImin2,xImax2
5631 integer,
parameter :: nsample=5
5632 integer :: iigrid,igrid,ir,it,ip
5633 integer :: ixI^L,ixO^L
5634 double precision,
allocatable :: rface(:),thetaface(:),phiface(:)
5635 double precision :: local_min1,local_max1,local_min2,local_max2
5636 double precision :: sph(1:3),xcent(1:2),wr,wt,wp
5638 local_min1=huge(one)
5639 local_max1=-huge(one)
5640 local_min2=huge(one)
5641 local_max2=-huge(one)
5643 do iigrid=1,igridstail
5644 igrid=igrids(iigrid)
5645 ^d&ixomin^d=ixmlo^d\
5646 ^d&ixomax^d=ixmhi^d\
5647 ^d&iximin^d=ixglo^d\
5648 ^d&iximax^d=ixghi^d\
5651 rface,thetaface,phiface)
5653 wr=dble(ir)/dble(nsample-1)
5654 sph(1)=(one-wr)*rface(ixomin1)+wr*rface(ixomax1+1)
5656 wt=dble(it)/dble(nsample-1)
5657 sph(2)=(one-wt)*thetaface(ixomin2)+wt*thetaface(ixomax2+1)
5659 wp=dble(ip)/dble(nsample-1)
5660 if (ir/=0 .and. ir/=nsample-1 .and. it/=0 .and. it/=nsample-1 .and. &
5661 ip/=0 .and. ip/=nsample-1) cycle
5662 sph(3)=(one-wp)*phiface(ixomin3)+wp*phiface(ixomax3+1)
5664 local_min1=min(local_min1,xcent(1))
5665 local_max1=max(local_max1,xcent(1))
5666 local_min2=min(local_min2,xcent(2))
5667 local_max2=max(local_max2,xcent(2))
5671 deallocate(rface,thetaface,phiface)
5674 call mpi_allreduce(local_min1,ximin1,1,mpi_double_precision,mpi_min,icomm,ierrmpi)
5675 call mpi_allreduce(local_max1,ximax1,1,mpi_double_precision,mpi_max,icomm,ierrmpi)
5676 call mpi_allreduce(local_min2,ximin2,1,mpi_double_precision,mpi_min,icomm,ierrmpi)
5677 call mpi_allreduce(local_max2,ximax2,1,mpi_double_precision,mpi_max,icomm,ierrmpi)
5678 if (ximin1>0.5d0*huge(one) .or. ximax1<-0.5d0*huge(one) .or. &
5679 ximin2>0.5d0*huge(one) .or. ximax2<-0.5d0*huge(one))
then
5680 call mpistop(
"sph_intersection could not determine image bounds")
5685 double precision,
intent(out) :: dxI
5687 select case(trim(dat_resolution_mode))
5693 call mpistop(
"unknown dat_resolution_mode")
5698 double precision,
intent(out) :: dxI
5700 double precision :: refine_factor,dr,dtheta,dphi,rmin,sin_theta_min
5702 refine_factor=dble(2**(refine_max_level-1))
5704 dxi=min(abs(xprobmax1-xprobmin1)/(dble(domain_nx1)*refine_factor),&
5705 abs(xprobmax2-xprobmin2)/(dble(domain_nx2)*refine_factor),&
5706 abs(xprobmax3-xprobmin3)/(dble(domain_nx3)*refine_factor))
5708 rmin=max(smalldouble,min(xprobmin1,xprobmax1))
5709 sin_theta_min=max(smalldouble,min(abs(sin(xprobmin2)),&
5710 abs(sin(xprobmax2))))
5711 dr=abs(xprobmax1-xprobmin1)/(dble(domain_nx1)*refine_factor)
5712 dtheta=abs(xprobmax2-xprobmin2)/(dble(domain_nx2)*refine_factor)
5713 dphi=abs(xprobmax3-xprobmin3)/(dble(domain_nx3)*refine_factor)
5714 dxi=min(dr,rmin*dtheta,rmin*sin_theta_min*dphi)
5716 call mpistop(
"nominal dat resolution needs Cartesian or spherical coordinates")
5719 if (dxi<=zero .or. dxi>half*huge(one))
then
5720 call mpistop(
"could not determine nominal dat-resolution image spacing")
5725 double precision,
intent(out) :: dxI
5727 integer :: iigrid,igrid,ixI^L,ixO^L,ix^D
5728 double precision :: local_min,global_min,dr,ds_theta,ds_phi,rval,theta
5731 do iigrid=1,igridstail
5732 igrid=igrids(iigrid)
5733 ^d&ixomin^d=ixmlo^d\
5734 ^d&ixomax^d=ixmhi^d\
5735 ^d&iximin^d=ixglo^d\
5736 ^d&iximax^d=ixghi^d\
5738 do ix1=ixomin1,ixomax1
5739 do ix2=ixomin2,ixomax2
5740 do ix3=ixomin3,ixomax3
5742 local_min=min(local_min,ps(igrid)%dx(ix^d,1),&
5743 ps(igrid)%dx(ix^d,2),ps(igrid)%dx(ix^d,3))
5745 rval=max(smalldouble,ps(igrid)%x(ix^d,1))
5746 theta=ps(igrid)%x(ix^d,2)
5747 dr=ps(igrid)%dx(ix^d,1)
5748 ds_theta=rval*ps(igrid)%dx(ix^d,2)
5749 ds_phi=rval*max(smalldouble,sin(theta))*ps(igrid)%dx(ix^d,3)
5750 local_min=min(local_min,dr,ds_theta,ds_phi)
5752 call mpistop(
"minimum dat resolution needs Cartesian or spherical coordinates")
5759 call mpi_allreduce(local_min,global_min,1,mpi_double_precision,mpi_min,&
5761 if (global_min<=zero .or. global_min>half*huge(one))
then
5762 call mpistop(
"could not determine minimum dat-resolution image spacing")
5773 integer,
intent(in) :: qunit
5775 character(20),
intent(in) :: datatype
5777 integer :: ix^D,numXI1,numXI2,numWI
5778 double precision :: xImin1,xImax1,xImin2,xImax2,xIcent1,xIcent2,dxI
5779 double precision,
allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:)
5780 double precision,
allocatable :: wI(:,:,:),wIs(:,:,:),EM(:,:),Dpl(:,:),Tau(:,:),EMthin(:,:),WLB(:,:,:)
5781 double precision :: vec_temp1(1:3),vec_temp2(1:3)
5782 double precision :: vec_z(1:3),vec_cor(1:3),xI_cor(1:2)
5783 double precision :: res,LOS_psi,r_max,r_loc
5786 character (30) :: ion
5787 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
5788 double precision :: arcsec,RHESSI_rsl,LASCO_rsl,pixel,R_occult,smallflux
5789 integer :: iigrid,igrid,i,j,numSI,iw
5790 logical :: emit,ray_image_global,has_thick_output
5804 ximin1=-abs(xprobmax1)
5805 ximin2=-abs(xprobmax1)
5806 ximax1=abs(xprobmax1)
5807 ximax2=abs(xprobmax1)
5812 if (ix1==1) vec_cor(1)=xprobmin1
5813 if (ix1==2) vec_cor(1)=xprobmax1
5815 if (ix2==1) vec_cor(2)=xprobmin2
5816 if (ix2==2) vec_cor(2)=xprobmax2
5818 if (ix3==1) vec_cor(3)=xprobmin3
5819 if (ix3==2) vec_cor(3)=xprobmax3
5822 r_loc=r_loc+(vec_cor(2)-
x_origin(2))**2
5823 r_loc=r_loc+(vec_cor(3)-
x_origin(3))**2
5825 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
5828 r_max=max(r_max,r_loc)
5832 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
5838 ximin1=min(ximin1,xi_cor(1))
5839 ximax1=max(ximax1,xi_cor(1))
5840 ximin2=min(ximin2,xi_cor(2))
5841 ximax2=max(ximax2,xi_cor(2))
5854 xicent1=(ximin1+ximax1)/2.d0
5855 xicent2=(ximin2+ximax2)/2.d0
5863 if (datatype==
'image_euv')
then
5867 else if (datatype==
'image_sxr')
then
5869 dxi=rhessi_rsl*arcsec
5871 else if (datatype==
'image_whitelight')
then
5882 call mpistop(
'Whitelight synthesis: instrument is not supported!')
5884 dxi=lasco_rsl*arcsec
5889 numxi1=8*ceiling((ximax1-xicent1)/dxi/8.d0)
5890 ximin1=xicent1-numxi1*dxi
5891 ximax1=xicent1+numxi1*dxi
5893 numxi2=8*ceiling((ximax2-xicent2)/dxi/8.d0)
5894 ximin2=xicent2-numxi2*dxi
5895 ximax2=xicent2+numxi2*dxi
5897 allocate(xi1(numxi1),xi2(numxi2),dxi1(numxi1),dxi2(numxi2))
5899 xi1(ix1)=ximin1+dxi*(ix1-
half)
5903 xi2(ix2)=ximin2+dxi*(ix2-
half)
5908 if (datatype==
'image_euv' .or. datatype==
'image_sxr')
then
5909 has_thick_output=datatype==
'image_euv' .and. trim(
radiation_transfer)==
'thick' .and. &
5912 if (datatype==
'image_euv')
then
5917 allocate(wi(numxi1,numxi2,numwi),wis(numxi1,numxi2,numwi),em(numxi1,numxi2))
5921 ray_image_global=.false.
5922 if (has_thick_output)
then
5923 allocate(tau(numxi1,numxi2),emthin(numxi1,numxi2))
5927 if (
slab .and. datatype==
'image_euv' .and. &
5929 ray_image_global=.true.
5930 allocate(dpl(numxi1,numxi2))
5939 do iigrid=1,igridstail; igrid=igrids(iigrid);
5942 else if (trim(
ray_method_active) ==
'spherical' .and. datatype ==
'image_euv')
then
5944 ray_image_global=.true.
5950 do iigrid=1,igridstail; igrid=igrids(iigrid);
5954 if (ray_image_global)
then
5955 if (has_thick_output)
then
5957 has_thick_output,tau=tau,euvthin=emthin,&
5958 cap_absorption=.true.)
5965 if (em(ix1,ix2)>smallflux) wis(ix1,ix2,1)=em(ix1,ix2)
5969 if (.not. ray_image_global)
then
5970 numsi=numxi1*numxi2*numwi
5971 call mpi_allreduce(wis,wi,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5979 call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
5980 if (
allocated(tau))
deallocate(tau)
5981 if (
allocated(emthin))
deallocate(emthin)
5982 deallocate(wi,wis,em)
5983 else if (datatype==
'image_whitelight')
then
5985 allocate(wi(numxi1,numxi2,numwi),wis(numxi1,numxi2,numwi),wlb(numxi1,numxi2,numwi))
5990 do iigrid=1,igridstail; igrid=igrids(iigrid);
5996 if (wlb(ix1,ix2,1)>smallflux)
then
5997 wis(ix1,ix2,1)=wlb(ix1,ix2,1)
5998 wis(ix1,ix2,2)=wlb(ix1,ix2,2)
6002 numsi=numxi1*numxi2*numwi
6003 call mpi_allreduce(wis,wi,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
6010 call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
6011 deallocate(wi,wis,wlb)
6014 deallocate(xi1,xi2,dxi1,dxi2)
6019 integer,
intent(in) :: igrid,numXI1,numXI2
6020 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
6021 double precision,
intent(in) :: dxI
6023 character(20),
intent(in) :: datatype
6024 double precision,
intent(inout) :: EM(numXI1,numXI2)
6026 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
6027 double precision :: xb^L,xd^D
6028 double precision,
allocatable :: flux(:^D&),opacity(:^D&)
6029 double precision :: res
6030 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
6031 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC
6032 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2)
6035 double precision :: logTe
6036 character (30) :: ion
6037 double precision :: lineCent
6038 double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
6039 double precision :: arcsec,pixel,RHESSI_rsl,area_1AU
6040 double precision :: aa,bb
6042 ^d&ixomin^d=ixmlo^d\
6043 ^d&ixomax^d=ixmhi^d\
6044 ^d&iximin^d=ixglo^d\
6045 ^d&iximax^d=ixghi^d\
6046 ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
6047 ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
6050 arcsec=7.25d5/unit_length
6052 arcsec=7.25d7/unit_length
6055 allocate(flux(ixi^s),opacity(ixi^s))
6056 if (datatype==
'image_euv')
then
6057 if (trim(emission_model)==
'pseudo_current')
then
6059 else if (trim(emission_model)==
'radio_ff')
then
6063 call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
6064 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
6066 call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
6067 pixel=spacersl*arcsec
6068 sigma0=sigma_psf*pixel
6069 else if (datatype==
'image_sxr')
then
6071 call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr)
6072 rhessi_rsl=2.3d0/instrument_resolution_factor
6074 pixel=rhessi_rsl*arcsec
6075 sigma0=sigma_psf*pixel
6080 {
do ix^d=ixomin^d,ixomax^d\}
6082 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi1(^d))/(dxi/2.d0)));
6083 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi2(^d))/(dxi/2.d0)));
6084 ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
6085 if (datatype==
'image_euv')
then
6087 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length*1.d2/dxi/dxi
6089 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length/dxi/dxi
6091 else if (datatype==
'image_sxr')
then
6093 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length**3/area_1au
6095 if (fluxsubc>smalldouble)
then
6097 {
do isubc^d=1,nsubc^d\}
6098 ^d&xsubc(^d)=ps(igrid)%x(ix^dd,^d)-half*ps(igrid)%dx(ix^dd,^d)+(isubc^d-half)*dxsubc^d;
6102 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
6103 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
6104 ixpmin1=max(1,ixp1-3)
6105 ixpmax1=min(ixp1+3,numxi1)
6106 ixpmin2=max(1,ixp2-3)
6107 ixpmax2=min(ixp2+3,numxi2)
6108 do ixp1=ixpmin1,ixpmax1
6109 do ixp2=ixpmin2,ixpmax2
6110 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6111 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6112 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6113 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6114 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
6115 em(ixp1,ixp2)=em(ixp1,ixp2)+fluxsubc*factor
6122 deallocate(flux,opacity)
6126 integer,
intent(in) :: igrid,numXI1,numXI2
6127 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
6128 double precision,
intent(in) :: dxI
6130 character(20),
intent(in) :: datatype
6131 double precision,
intent(inout) :: EM(numXI1,numXI2)
6133 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
6134 double precision,
allocatable :: flux(:^D&),Ne(:^D&),opacity(:^D&)
6135 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
6136 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC,RsubC
6137 double precision :: TBsubC,PBsubC
6138 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2),xSubC_car(1:3)
6139 double precision :: R_thick,dotp,dvolume,R_occult,Rc
6140 double precision :: dxl(1:3),x_sph(1:3),dx_sph(1:3)
6141 double precision :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
6142 logical :: sun_back_side,emit
6145 double precision :: logTe
6146 character (30) :: ion
6147 double precision :: lineCent
6148 double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
6149 double precision :: RHESSI_rsl,area_1AU,arcsec,pixel
6151 ^d&ixomin^d=ixmlo^d;
6152 ^d&ixomax^d=ixmhi^d;
6153 ^d&iximin^d=ixglo^d;
6154 ^d&iximax^d=ixghi^d;
6157 arcsec=7.25d5/unit_length
6159 arcsec=7.25d7/unit_length
6162 allocate(flux(ixi^s),opacity(ixi^s))
6163 if (datatype==
'image_euv')
then
6164 if (trim(emission_model)==
'pseudo_current')
then
6166 else if (trim(emission_model)==
'radio_ff')
then
6170 call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
6171 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
6173 call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
6174 pixel=spacersl*arcsec
6175 sigma0=sigma_psf*pixel
6176 else if (datatype==
'image_sxr')
then
6178 call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr)
6179 rhessi_rsl=2.3d0/instrument_resolution_factor
6181 pixel=rhessi_rsl*arcsec
6182 sigma0=sigma_psf*pixel
6187 r_thick=r_opt_thick*const_rsun/unit_length
6188 {
do ix^d=ixomin^d,ixomax^d\}
6189 x_sph(1:3)=ps(igrid)%x(ix^d,1:3)
6190 dx_sph(1:3)=ps(igrid)%dx(ix^d,1:3)
6192 dxl(2)=x_sph(1)*dx_sph(2)
6193 dxl(3)=x_sph(1)*dsin(x_sph(2))*dx_sph(3)
6198 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
6200 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
6202 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
6204 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
6206 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
6208 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
6213 xsubc(1)=x_sph(1)-half*dx_sph(1)+(isubc1-half)*dx_sph(1)/nsubc1
6215 dxsubc1=dx_sph(1)/nsubc1
6218 xsubc(2)=x_sph(2)-half*dx_sph(2)+(isubc2-half)*dx_sph(2)/nsubc2
6219 dxsubc2=xsubc(1)*dx_sph(2)/nsubc2
6220 dxsubc3=xsubc(1)*dsin(xsubc(2))*dx_sph(3)/nsubc3
6221 dvolume=dxsubc1*dxsubc2*dxsubc3
6222 if (datatype==
'image_euv')
then
6224 fluxsubc=flux(ix^d)*dvolume*unit_length*1.d2/dxi/dxi
6226 fluxsubc=flux(ix^d)*dvolume*unit_length/dxi/dxi
6228 else if (datatype==
'image_sxr')
then
6230 fluxsubc=flux(ix^d)*dvolume*unit_length**3/area_1au
6233 if (fluxsubc>smalldouble)
then
6236 xsubc(3)=x_sph(3)-half*dx_sph(3)+(isubc3-half)*dx_sph(3)/nsubc3
6238 rc=dsqrt(xcent(1)**2+xcent(2)**2)
6243 sun_back_side=.true.
6244 if (dotp<0.d0) sun_back_side=.false.
6246 if (sun_back_side)
then
6248 if (rc>r_thick) emit=.true.
6251 if (xsubc(1)<=r_thick) emit=.false.
6257 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
6258 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
6259 ixpmin1=max(1,ixp1-3)
6260 ixpmax1=min(ixp1+3,numxi1)
6261 ixpmin2=max(1,ixp2-3)
6262 ixpmax2=min(ixp2+3,numxi2)
6263 do ixp1=ixpmin1,ixpmax1
6264 do ixp2=ixpmin2,ixpmax2
6265 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6266 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6267 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6268 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6269 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
6270 em(ixp1,ixp2)=em(ixp1,ixp2)+fluxsubc*factor
6280 deallocate(flux,opacity)
6287 integer,
intent(in) :: igrid,numXI1,numXI2,numWI
6288 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
6289 double precision,
intent(in) :: dxI
6291 character(20),
intent(in) :: datatype
6292 double precision,
intent(inout) :: WLB(numXI1,numXI2,numWI)
6294 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
6295 double precision,
allocatable :: flux(:^D&),Ne(:^D&),nH_dummy(:^D&)
6296 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
6297 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC,RsubC
6298 double precision :: sigma_PSF,sigma0,arcsec,pixel,LASCO_rsl
6299 double precision :: A,B,C,D,Rc,Ne0,TBsubC,PBsubC,factor
6300 double precision :: R_thick,dotp,dvolume,R_occult
6301 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2),xSubC_car(1:3)
6302 double precision :: dxl(1:3),x_sph(1:3),dx_sph(1:3)
6303 double precision :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
6306 ^d&ixomin^d=ixmlo^d;
6307 ^d&ixomax^d=ixmhi^d;
6308 ^d&iximin^d=ixglo^d;
6309 ^d&iximax^d=ixghi^d;
6312 arcsec=7.25d5/unit_length
6314 arcsec=7.25d7/unit_length
6318 allocate(nh_dummy(ixi^s))
6319 if (whitelight_instrument==
'LASCO/C1')
then
6320 lasco_rsl=5.6d0/instrument_resolution_factor
6322 else if (whitelight_instrument==
'LASCO/C2')
then
6323 lasco_rsl=11.4d0/instrument_resolution_factor
6325 else if (whitelight_instrument==
'LASCO/C3')
then
6326 lasco_rsl=56.d0/instrument_resolution_factor
6329 if (r_occultor>1.d0) r_occult=r_occultor
6330 r_occult=r_occult*const_rsun/unit_length
6331 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,ne)
6333 call eos%get_ne_nH(ixi^l, ixo^l, ps(igrid)%w, ps(igrid)%x, ne, nh_dummy)
6335 pixel=lasco_rsl*arcsec
6336 sigma0=sigma_psf*pixel
6339 r_thick=r_opt_thick*const_rsun/unit_length
6340 {
do ix^d=ixomin^d,ixomax^d\}
6341 x_sph(1:3)=ps(igrid)%x(ix^d,1:3)
6342 dx_sph(1:3)=ps(igrid)%dx(ix^d,1:3)
6344 dxl(2)=x_sph(1)*dx_sph(2)
6345 dxl(3)=x_sph(1)*dsin(x_sph(2))*dx_sph(3)
6346 ne0=ne(ix^d)*unit_numberdensity
6351 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
6353 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
6355 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
6357 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
6359 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
6361 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
6366 xsubc(1)=x_sph(1)-half*dx_sph(1)+(isubc1-half)*dx_sph(1)/nsubc1
6368 dxsubc1=dx_sph(1)/nsubc1
6372 xsubc(2)=x_sph(2)-half*dx_sph(2)+(isubc2-half)*dx_sph(2)/nsubc2
6373 dxsubc2=xsubc(1)*dx_sph(2)/nsubc2
6374 dxsubc3=xsubc(1)*dsin(xsubc(2))*dx_sph(3)/nsubc3
6375 dvolume=dxsubc1*dxsubc2*dxsubc3
6378 xsubc(3)=x_sph(3)-half*dx_sph(3)+(isubc3-half)*dx_sph(3)/nsubc3
6380 rc=dsqrt(xcent(1)**2+xcent(2)**2)
6383 if (rc>r_occult)
then
6387 tbsubc=tbsubc*dvolume*unit_length/dxi/dxi
6388 pbsubc=pbsubc*dvolume*unit_length/dxi/dxi
6389 if (tbsubc<1.d-20) emit=.false.
6394 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
6395 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
6396 ixpmin1=max(1,ixp1-3)
6397 ixpmax1=min(ixp1+3,numxi1)
6398 ixpmin2=max(1,ixp2-3)
6399 ixpmax2=min(ixp2+3,numxi2)
6400 do ixp1=ixpmin1,ixpmax1
6401 do ixp2=ixpmin2,ixpmax2
6402 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6403 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6404 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6405 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6406 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
6407 wlb(ixp1,ixp2,1)=wlb(ixp1,ixp2,1)+tbsubc*factor
6408 wlb(ixp1,ixp2,2)=wlb(ixp1,ixp2,2)+pbsubc*factor
6418 deallocate(nh_dummy)
6425 double precision,
intent(in) :: Rl
6426 double precision,
intent(inout) :: A,B,C,D
6428 double precision :: sinO,cosO,sinO2,cosO2,tmp
6433 coso=abs(dsqrt(coso2))
6434 tmp=log((1.d0+sino)/coso)
6436 b=-(1.d0-3.d0*sino2-(coso2/sino)*(1.d0+3.d0*sino2)*tmp)/8.d0
6437 c=4.d0/3.d0-coso-coso*coso2/3.d0
6438 d=(5.d0+sino2-(coso2/sino)*(5.d0-sino2)*tmp)/8.d0
6444 double precision,
intent(in) :: Rl,Rin,Ne,A,B,C,D
6445 double precision,
intent(inout) :: fluxTB,fluxPB
6447 double precision :: const,u,Bt,Br,PB,TB,sinchi2
6450 const=1.24878d-25/(1.d0-u/3.d0)
6452 bt=const*(c+u*(d-c))
6453 pb=const*sinchi2*((a+u*(b-a)))
6462 double precision,
intent(in) :: x_sph(1:3)
6463 double precision,
intent(inout) :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
6465 unitv_r(1)=dsin(x_sph(2))*dcos(x_sph(3))
6466 unitv_r(2)=dsin(x_sph(2))*dsin(x_sph(3))
6467 unitv_r(3)=dcos(x_sph(2))
6468 unitv_theta(1)=dcos(x_sph(2))*dcos(x_sph(3))
6469 unitv_theta(2)=dcos(x_sph(2))*dsin(x_sph(3))
6470 unitv_theta(3)=-dsin(x_sph(2))
6471 unitv_phi(1)=-dsin(x_sph(3))
6472 unitv_phi(2)=dcos(x_sph(3))
6477 subroutine output_data(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,datatype)
6481 integer,
intent(in) :: qunit,nXO1,nXO2,nWO
6482 double precision,
intent(in) :: dxO1(nxO1),dxO2(nxO2)
6483 double precision,
intent(in) :: xO1(nXO1),xO2(nxO2)
6484 double precision,
intent(inout) :: wO(nXO1,nXO2,nWO)
6485 character(20),
intent(in) :: datatype
6487 integer :: nPiece,nP1,nP2,nC1,nC2,nWC
6488 integer :: piece_nmax1,piece_nmax2,ix1,ix2,j,ipc,ixc1,ixc2
6489 double precision :: uniform_tol
6490 double precision,
allocatable :: xC(:,:,:,:),wC(:,:,:,:),dxC(:,:,:,:)
6497 if (abs(wo(ix1,ix2,j))<smalldouble) wo(ix1,ix2,j)=zero
6504 if (datatype==
'image_euv' .or. datatype==
'image_sxr')
then
6506 piece_nmax1=block_nx2
6507 piece_nmax2=block_nx3
6509 piece_nmax1=block_nx3
6510 piece_nmax2=block_nx1
6512 piece_nmax1=block_nx1
6513 piece_nmax2=block_nx2
6515 else if (datatype==
'spectrum_euv')
then
6518 piece_nmax2=block_nx1
6520 piece_nmax2=block_nx2
6522 piece_nmax2=block_nx3
6529 loopn1:
do j=piece_nmax1,1,-1
6530 if(mod(nxo1,j)==0)
then
6535 loopn2:
do j=piece_nmax2,1,-1
6536 if(mod(nxo2,j)==0)
then
6549 case(
'EIvtuCCmpi',
'ESvtuCCmpi',
'SIvtuCCmpi',
'WIvtuCCmpi')
6551 allocate(xc(npiece,nc1,nc2,2))
6552 allocate(dxc(npiece,nc1,nc2,2))
6553 allocate(wc(npiece,nc1,nc2,nwo))
6557 ix1=mod(ipc-1,np1)*nc1+ixc1
6558 ix2=floor(1.0*(ipc-1)/np1)*nc2+ixc2
6559 xc(ipc,ixc1,ixc2,1)=xo1(ix1)
6560 xc(ipc,ixc1,ixc2,2)=xo2(ix2)
6561 dxc(ipc,ixc1,ixc2,1)=dxo1(ix1)
6562 dxc(ipc,ixc1,ixc2,2)=dxo2(ix2)
6564 wc(ipc,ixc1,ixc2,j)=wo(ix1,ix2,j)
6571 deallocate(xc,dxc,wc)
6572 case(
'EIvtiCCmpi',
'ESvtiCCmpi',
'SIvtiCCmpi',
'WIvtiCCmpi')
6574 (maxval(abs(dxo1(:)-dxo1(1)))>uniform_tol*max(one,abs(dxo1(1))) .or. &
6575 maxval(abs(dxo2(:)-dxo2(1)))>uniform_tol*max(one,abs(dxo2(1)))))
then
6576 call mpistop(
"vti needs uniform dat-resolution image grids")
6578 call write_image_vticc(qunit,xo1,xo2,dxo1,dxo2,wo,nxo1,nxo2,nwo,nc1,nc2)
6581 call mpistop(
"Error in synthesize emission: Unknown convert_type")
6587 subroutine write_image_vticc(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,nC1,nC2)
6591 integer,
intent(in) :: qunit,nXO1,nXO2,nWO,nC1,nC2
6592 double precision,
intent(in) :: xO1(nXO1),xO2(nxO2)
6593 double precision,
intent(in) :: dxO1(nxO1),dxO2(nxO2)
6594 double precision,
intent(in) :: wO(nXO1,nXO2,nWO)
6596 double precision :: origin(1:3), spacing(1:3)
6597 integer :: wholeExtent(1:6)
6599 integer :: ixC1,ixC2
6603 character (70) :: subname,wname,vname,nameL,nameS
6604 character (len=std_len) :: filename
6605 logical :: sph_datres_no_doppler
6608 origin(1)=xo1(1)-0.5d0*dxo1(1)
6609 origin(2)=xo2(1)-0.5d0*dxo2(1)
6620 inquire(qunit,opened=fileopen)
6621 if(.not.fileopen)
then
6626 write(filename,
'(a,i4.4,a)') trim(
filename_euv),filenr,
".vti"
6628 write(filename,
'(a,i4.4,a)') trim(
filename_sxr),filenr,
".vti"
6634 open(qunit,file=filename,status=
'unknown',form=
'formatted')
6638 write(qunit,
'(a)')
'<?xml version="1.0"?>'
6639 write(qunit,
'(a)',advance=
'no')
'<VTKFile type="ImageData"'
6640 write(qunit,
'(a)')
' version="0.1" byte_order="LittleEndian">'
6641 write(qunit,
'(a,3(1pe14.6),a,6(i10),a,3(1pe14.6),a)')
' <ImageData Origin="',&
6642 origin,
'" WholeExtent="',wholeextent,
'" Spacing="',spacing,
'">'
6644 write(qunit,
'(a)')
'<FieldData>'
6645 write(qunit,
'(2a)')
'<DataArray type="Float32" Name="TIME" ',&
6646 'NumberOfTuples="1" format="ascii">'
6648 write(qunit,
'(a)')
'</DataArray>'
6649 write(qunit,
'(a)')
'</FieldData>'
6651 write(qunit,
'(a,6(i10),a)')
'<Piece Extent="',wholeextent,
'">'
6652 write(qunit,
'(a)')
'<CellData>'
6663 if (trim(
emission_model)==
'pseudo_current' .and. iw==1) vname=
'pseudo_current'
6664 if (trim(
emission_model)==
'radio_ff' .and. iw==1) vname=
'radio_brightness_temperature'
6666 if (iw==2 .and.
dat_resolution .and. (.not. sph_datres_no_doppler) .and. &
6672 ((
dat_resolution .and. ((sph_datres_no_doppler .and. iw==2) .or. &
6673 ((.not. sph_datres_no_doppler) .and. iw==3))) .or. &
6687 vname=
'absorption_fraction'
6698 if (iw==1)
write(vname,
'(a)')
'B'
6699 if (iw==2)
write(vname,
'(a)')
'pB'
6707 write(qunit,
'(a,a,a)')&
6708 '<DataArray type="Float64" Name="',trim(vname),
'" format="ascii">'
6709 write(qunit,
'(200(1pe14.6))') ((wo(ixc1,ixc2,iw),ixc1=1,nxo1),ixc2=1,nxo2)
6710 write(qunit,
'(a)')
'</DataArray>'
6712 write(qunit,
'(a)')
'</CellData>'
6713 write(qunit,
'(a)')
'</Piece>'
6715 write(qunit,
'(a)')
'</ImageData>'
6716 write(qunit,
'(a)')
'</VTKFile>'
6726 integer,
intent(in) :: qunit
6727 integer,
intent(in) :: nPiece,nC1,nC2,nWC
6728 double precision,
intent(in) :: xC(nPiece,nC1,nC2,2),dxC(nPiece,nc1,nc2,2)
6729 double precision,
intent(in) :: wC(nPiece,nC1,nC2,nWC)
6730 character(20),
intent(in) :: datatype
6733 double precision :: xP(nPiece,nC1+1,nC2+1,2)
6736 character (70) :: subname,wname,vname,nameL,nameS
6737 character (len=std_len) :: filename
6738 integer :: ixC1,ixC2,ixP,ix1,ix2,j
6739 integer :: nc,np,icel,VTK_type
6740 logical :: sph_datres_no_doppler
6751 if (ix1<np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1,1,1)-0.5d0*dxc(ixp,ix1,1,1)
6752 if (ix1==np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1-1,1,1)+0.5d0*dxc(ixp,ix1-1,1,1)
6753 if (ix2<np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2,2)-0.5d0*dxc(ixp,1,ix2,2)
6754 if (ix2==np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2-1,2)+0.5d0*dxc(ixp,1,ix2-1,2)
6759 inquire(qunit,opened=fileopen)
6760 if(.not.fileopen)
then
6764 if (datatype==
'image_euv')
then
6765 write(filename,
'(a,i4.4,a)') trim(
filename_euv),filenr,
".vtu"
6766 else if (datatype==
'image_sxr')
then
6767 write(filename,
'(a,i4.4,a)') trim(
filename_sxr),filenr,
".vtu"
6768 else if (datatype==
'image_whitelight')
then
6770 else if (datatype==
'spectrum_euv')
then
6773 open(qunit,file=filename,status=
'unknown',form=
'formatted')
6776 write(qunit,
'(a)')
'<?xml version="1.0"?>'
6777 write(qunit,
'(a)',advance=
'no')
'<VTKFile type="UnstructuredGrid"'
6778 write(qunit,
'(a)')
' version="0.1" byte_order="LittleEndian">'
6779 write(qunit,
'(a)')
'<UnstructuredGrid>'
6780 write(qunit,
'(a)')
'<FieldData>'
6781 write(qunit,
'(2a)')
'<DataArray type="Float32" Name="TIME" ',&
6782 'NumberOfTuples="1" format="ascii">'
6784 write(qunit,
'(a)')
'</DataArray>'
6785 write(qunit,
'(a)')
'</FieldData>'
6787 write(qunit,
'(a,i7,a,i7,a)') &
6788 '<Piece NumberOfPoints="',np,
'" NumberOfCells="',nc,
'">'
6789 write(qunit,
'(a)')
'<CellData>'
6791 if (datatype==
'image_euv')
then
6800 if (trim(
emission_model)==
'pseudo_current') vname=
'pseudo_current'
6801 if (trim(
emission_model)==
'radio_ff') vname=
'radio_brightness_temperature'
6804 if (j==2 .and.
dat_resolution .and. (.not. sph_datres_no_doppler) .and. &
6810 ((
dat_resolution .and. ((sph_datres_no_doppler .and. j==2) .or. &
6811 ((.not. sph_datres_no_doppler) .and. j==3))) .or. &
6825 vname=
'absorption_fraction'
6827 else if (datatype==
'image_sxr')
then
6835 else if (datatype==
'image_whitelight')
then
6836 write(vname,
'(a)')
'whitelight'
6837 else if (datatype==
'spectrum_euv')
then
6844 write(qunit,
'(a,a,a)')&
6845 '<DataArray type="Float64" Name="',trim(vname),
'" format="ascii">'
6846 write(qunit,
'(200(1pe14.6))') ((wc(ixp,ixc1,ixc2,j),ixc1=1,nc1),ixc2=1,nc2)
6847 write(qunit,
'(a)')
'</DataArray>'
6849 write(qunit,
'(a)')
'</CellData>'
6850 write(qunit,
'(a)')
'<Points>'
6851 write(qunit,
'(a)')
'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
6856 write(qunit,
'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
6858 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
6860 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
6864 write(qunit,
'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
6866 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
6868 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
6871 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
6875 write(qunit,
'(a)')
'</DataArray>'
6876 write(qunit,
'(a)')
'</Points>'
6878 write(qunit,
'(a)')
'<Cells>'
6879 write(qunit,
'(a)')
'<DataArray type="Int32" Name="connectivity" format="ascii">'
6882 write(qunit,
'(4(i7))') ix1-1+(ix2-1)*np1,ix1+(ix2-1)*np1,&
6883 ix1-1+ix2*np1,ix1+ix2*np1
6886 write(qunit,
'(a)')
'</DataArray>'
6888 write(qunit,
'(a)')
'<DataArray type="Int32" Name="offsets" format="ascii">'
6890 write(qunit,
'(i7)') icel*(2**2)
6892 write(qunit,
'(a)')
'</DataArray>'
6894 write(qunit,
'(a)')
'<DataArray type="Int32" Name="types" format="ascii">'
6898 write(qunit,
'(i2)') vtk_type
6900 write(qunit,
'(a)')
'</DataArray>'
6901 write(qunit,
'(a)')
'</Cells>'
6902 write(qunit,
'(a)')
'</Piece>'
6904 write(qunit,
'(a)')
'</UnstructuredGrid>'
6905 write(qunit,
'(a)')
'</VTKFile>'
6911 double precision,
intent(in) :: vec1(1:3),vec2(1:3)
6912 double precision,
intent(out) :: res
6914 res=vec1(1)*vec2(1)+vec1(2)*vec2(2)+vec1(3)*vec2(3)
6919 double precision,
intent(in) :: vec_in1(1:3),vec_in2(1:3)
6920 double precision,
intent(out) :: vec_out(1:3)
6922 vec_out(1)=vec_in1(2)*vec_in2(3)-vec_in1(3)*vec_in2(2)
6923 vec_out(2)=vec_in1(3)*vec_in2(1)-vec_in1(1)*vec_in2(3)
6924 vec_out(3)=vec_in1(1)*vec_in2(2)-vec_in1(2)*vec_in2(1)
6930 double precision :: LOS_psi
6931 double precision :: vec_car(1:3),vec_z(1:3),vec_temp1(1:3),vec_temp2(1:3)
6932 double precision :: vec_LOS_sph(1:3),vec_xI1_sph(1:3),vec_xI2_sph(1:3)
6936 vec_los(2)=dpi*los_theta/180.d0
6947 if (los_theta==zero)
then
6949 vec_temp1(2)=dpi/2.d0
6950 vec_temp1(3)=dpi*los_phi/180.d0
6964 los_psi=dpi*image_rotate/180.d0
6965 vec_xi1=vec_temp1*cos(los_psi)-vec_temp2*sin(los_psi)
6966 vec_xi2=vec_temp2*cos(los_psi)+vec_temp1*sin(los_psi)
6977 vec_los_sph(2:3)=vec_los_sph(2:3)*180.d0/dpi
6978 vec_xi1_sph(2:3)=vec_xi1_sph(2:3)*180.d0/dpi
6979 vec_xi2_sph(2:3)=vec_xi2_sph(2:3)*180.d0/dpi
6981 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),
']'
6982 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),
']'
6983 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),
']'
6989 double precision,
intent(in) :: vec_sph(1:3)
6990 double precision,
intent(inout) :: vec_car(1:3)
6992 vec_car(1)=vec_sph(1)*dsin(vec_sph(2))*dcos(vec_sph(3))
6993 vec_car(2)=vec_sph(1)*dsin(vec_sph(2))*dsin(vec_sph(3))
6994 vec_car(3)=vec_sph(1)*dcos(vec_sph(2))
7000 double precision,
intent(in) :: vec_car(1:3)
7001 double precision,
intent(inout) :: vec_sph(1:3)
7003 vec_sph(1)=dsqrt(vec_car(1)**2+vec_car(2)**2+vec_car(3)**2)
7004 vec_sph(2)=dacos(vec_car(3)/vec_sph(1))
7005 vec_sph(3)=atan2(vec_car(2),vec_car(1))
7011 double precision :: LOS_psi
7012 double precision :: vec_z(1:3),vec_temp1(1:3),vec_temp2(1:3)
7015 vec_los(1)=-cos(dpi*los_phi/180.d0)*sin(dpi*los_theta/180.d0)
7016 vec_los(2)=-sin(dpi*los_phi/180.d0)*sin(dpi*los_theta/180.d0)
7017 vec_los(3)=-cos(dpi*los_theta/180.d0)
7023 if (los_theta==zero)
then
7024 vec_xi1(1)=cos(dpi*los_phi/180.d0)
7025 vec_xi1(2)=sin(dpi*los_phi/180.d0)
7033 los_psi=dpi*image_rotate/180.d0
7034 vec_xi1=vec_temp1*cos(los_psi)-vec_temp2*sin(los_psi)
7035 vec_xi2=vec_temp2*cos(los_psi)+vec_temp1*sin(los_psi)
7049 double precision,
intent(in) :: x_3D_sph(1:3)
7050 double precision,
intent(inout) :: x_image(1:2)
7051 double precision :: res,res_origin
7052 double precision :: x_3D(1:3)
7063 double precision,
intent(in) :: x_3D(1:3)
7064 double precision,
intent(inout) :: x_image(1:2)
7065 double precision :: res,res_origin
7069 x_image(1)=res-res_origin
7072 x_image(2)=res-res_origin
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
double precision, parameter const_rsun
double precision, parameter kb_cgs
Boltzmann constant in cgs.
double precision, parameter half
double precision, parameter one
double precision, parameter dpi
Pi.
double precision, parameter zero
some frequently used numbers
double precision, parameter smalldouble
double precision, parameter mp_cgs
Proton mass in cgs.
double precision, parameter const_c
universal constants as specified in cgs units
PI (partial-ionisation) ionisation-degree backend for the eos% family.
subroutine, public ionization_state_tp(t, p, rfactor, iz_h, iz_he)
Equation of state for AMRVAC, handled through a single eos_container object.
Module with geometry-related routines (e.g., divergence, curl).
integer, parameter spherical
integer, parameter cartesian
subroutine curlvector(qvec, ixil, ixol, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
This module contains definitions of global parameters and variables and some generic functions/subrou...
character(len=std_len) filename_sxr
Base file name for synthetic SXR emission output.
integer spectrum_wl
wave length for spectrum
integer ixghi
Upper index of grid block arrays.
logical activate_unit_arcsec
use arcsec as length unit of images/spectra
character(len=std_len) filename_spectrum
Base file name for synthetic EUV spectrum output.
double precision global_time
The global simulation time.
logical output_absorption_fraction
output absorption fraction for thick/thin EUV synthesis when available
double precision radio_beam_fwhm
Gaussian radio beam full width at half maximum in arcsec.
integer snapshotini
Resume from the snapshot with this index.
character(len=std_len) filename_euv
Base file name for synthetic EUV emission output.
logical instrument_postprocess
Post-process dat-resolution EUV images onto the instrument pixel grid.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) filename_whitelight
Base file name for synthetic white light.
character(len=std_len) convert_type
Which format to use when converting.
integer, parameter ndim
Number of spatial dimensions for grid variables.
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)