19 double precision ::
t_aia(1:101)
22 double precision ::
f_335(1:101)
35 data t_aia / 4. , 4.05, 4.1, 4.15, 4.2, 4.25, 4.3, 4.35, &
36 4.4, 4.45, 4.5, 4.55, 4.6, 4.65, 4.7, 4.75, &
37 4.8, 4.85, 4.9, 4.95, 5. , 5.05, 5.1, 5.15, &
38 5.2, 5.25, 5.3, 5.35, 5.4, 5.45, 5.5, 5.55, &
39 5.6, 5.65, 5.7, 5.75, 5.8, 5.85, 5.9, 5.95, &
40 6. , 6.05, 6.1, 6.15, 6.2, 6.25, 6.3, 6.35, &
41 6.4, 6.45, 6.5, 6.55, 6.6, 6.65, 6.7, 6.75, &
42 6.8, 6.85, 6.9, 6.95, 7. , 7.05, 7.1, 7.15, &
43 7.2, 7.25, 7.3, 7.35, 7.4, 7.45, 7.5, 7.55, &
44 7.6, 7.65, 7.7, 7.75, 7.8, 7.85, 7.9, 7.95, &
45 8. , 8.05, 8.1, 8.15, 8.2, 8.25, 8.3, 8.35, &
46 8.4, 8.45, 8.5, 8.55, 8.6, 8.65, 8.7, 8.75, &
47 8.8, 8.85, 8.9, 8.95, 9. /
49 data f_94 / 4.25022959
d-37, 4.35880298
d-36, 3.57054296
d-35, 2.18175426
d-34, &
50 8.97592571
d-34, 2.68512961
d-33, 7.49559346
d-33, 2.11603751
d-32, &
51 5.39752853
d-32, 1.02935904
d-31, 1.33822307
d-31, 1.40884290
d-31, &
52 1.54933156
d-31, 2.07543102
d-31, 3.42026227
d-31, 6.31171444
d-31, &
53 1.16559416
d-30, 1.95360497
d-30, 2.77818735
d-30, 3.43552578
d-30, &
54 4.04061803
d-30, 4.75470982
d-30, 5.65553769
d-30, 6.70595782
d-30, &
55 7.80680354
d-30, 8.93247715
d-30, 1.02618156
d-29, 1.25979030
d-29, &
56 1.88526483
d-29, 3.62448572
d-29, 7.50553279
d-29, 1.42337571
d-28, &
57 2.37912813
d-28, 3.55232305
d-28, 4.84985757
d-28, 6.20662827
d-28, &
58 7.66193687
d-28, 9.30403645
d-28, 1.10519802
d-27, 1.25786927
d-27, &
59 1.34362634
d-27, 1.33185242
d-27, 1.22302081
d-27, 1.05677973
d-27, &
60 9.23064720
d-28, 8.78570994
d-28, 8.02397416
d-28, 5.87681142
d-28, &
61 3.82272695
d-28, 3.11492649
d-28, 3.85736090
d-28, 5.98893519
d-28, &
62 9.57553548
d-28, 1.46650267
d-27, 2.10365847
d-27, 2.79406671
d-27, &
63 3.39420087
d-27, 3.71077520
d-27, 3.57296767
d-27, 2.95114380
d-27, &
64 2.02913103
d-27, 1.13361825
d-27, 5.13405629
d-28, 2.01305089
d-28, &
65 8.15781482
d-29, 4.28366817
d-29, 3.08701543
d-29, 2.68693906
d-29, &
66 2.51764203
d-29, 2.41773103
d-29, 2.33996083
d-29, 2.26997246
d-29, &
67 2.20316143
d-29, 2.13810001
d-29, 2.07424438
d-29, 2.01149189
d-29, &
68 1.94980213
d-29, 1.88917920
d-29, 1.82963583
d-29, 1.77116920
d-29, &
69 1.71374392
d-29, 1.65740593
d-29, 1.60214447
d-29, 1.54803205
d-29, &
70 1.49510777
d-29, 1.44346818
d-29, 1.39322305
d-29, 1.34441897
d-29, &
71 1.29713709
d-29, 1.25132618
d-29, 1.20686068
d-29, 1.14226584
d-29, &
72 1.09866413
d-29, 1.05635524
d-29, 1.01532444
d-29, 9.75577134
d-30, &
73 9.37102736
d-30, 8.99873335
d-30, 8.63860172
d-30, 8.29051944
d-30, &
76 data f_131 / 3.18403601
d-37, 3.22254703
d-36, 2.61657920
d-35, &
77 1.59575286
d-34, 6.65779556
d-34, 2.07015132
d-33, &
78 6.05768615
d-33, 1.76074833
d-32, 4.52633001
d-32, &
79 8.57121883
d-32, 1.09184271
d-31, 1.10207963
d-31, &
80 1.11371658
d-31, 1.29105226
d-31, 1.80385897
d-31, &
81 3.27295431
d-31, 8.92002136
d-31, 3.15214579
d-30, &
82 9.73440787
d-30, 2.22709702
d-29, 4.01788984
d-29, &
83 6.27471832
d-29, 8.91764995
d-29, 1.18725647
d-28, &
84 1.52888040
d-28, 2.05082946
d-28, 3.47651873
d-28, &
85 8.80482184
d-28, 2.66533063
d-27, 7.05805149
d-27, &
86 1.46072515
d-26, 2.45282476
d-26, 3.55303726
d-26, &
87 4.59075911
d-26, 5.36503515
d-26, 5.68444094
d-26, &
88 5.47222296
d-26, 4.81119761
d-26, 3.85959059
d-26, &
89 2.80383406
d-26, 1.83977650
d-26, 1.11182849
d-26, &
90 6.50748885
d-27, 3.96843481
d-27, 2.61876319
d-27, &
91 1.85525324
d-27, 1.39717024
d-27, 1.11504283
d-27, &
92 9.38169611
d-28, 8.24801234
d-28, 7.43331919
d-28, &
93 6.74537063
d-28, 6.14495760
d-28, 5.70805277
d-28, &
94 5.61219786
d-28, 6.31981777
d-28, 9.19747307
d-28, &
95 1.76795732
d-27, 3.77985446
d-27, 7.43166191
d-27, &
96 1.19785603
d-26, 1.48234676
d-26, 1.36673114
d-26, &
97 9.61047146
d-27, 5.61209353
d-27, 3.04779780
d-27, &
98 1.69378976
d-27, 1.02113491
d-27, 6.82223774
d-28, &
99 5.02099099
d-28, 3.99377760
d-28, 3.36279037
d-28, &
100 2.94767378
d-28, 2.65740865
d-28, 2.44396277
d-28, &
101 2.28003967
d-28, 2.14941419
d-28, 2.04178995
d-28, &
102 1.95031045
d-28, 1.87011994
d-28, 1.79777869
d-28, &
103 1.73093957
d-28, 1.66795789
d-28, 1.60785455
d-28, &
104 1.55002399
d-28, 1.49418229
d-28, 1.44022426
d-28, &
105 1.38807103
d-28, 1.33772767
d-28, 1.28908404
d-28, &
106 1.24196208
d-28, 1.17437501
d-28, 1.12854330
d-28, &
107 1.08410498
d-28, 1.04112003
d-28, 9.99529904
d-29, &
108 9.59358806
d-29, 9.20512291
d-29, 8.83009123
d-29, &
109 8.46817043
d-29, 8.11921928
d-29 /
111 data f_171 / 2.98015581
d-42, 1.24696230
d-40, 3.37614652
d-39, 5.64103034
d-38, &
112 5.20550266
d-37, 2.77785939
d-36, 1.16283616
d-35, 6.50007689
d-35, &
113 9.96177399
d-34, 1.89586076
d-32, 2.10982799
d-31, 1.36946479
d-30, &
114 6.27396553
d-30, 2.29955134
d-29, 7.13430211
d-29, 1.91024282
d-28, &
115 4.35358848
d-28, 7.94807808
d-28, 1.07431875
d-27, 1.08399488
d-27, &
116 9.16212938
d-28, 7.34715770
d-28, 6.59246382
d-28, 9.13541375
d-28, &
117 2.05939035
d-27, 5.08206555
d-27, 1.10148083
d-26, 2.01884662
d-26, &
118 3.13578384
d-26, 4.14367719
d-26, 5.36067711
d-26, 8.74170213
d-26, &
119 1.64161233
d-25, 2.94587860
d-25, 4.76298332
d-25, 6.91765639
d-25, &
120 9.08825111
d-25, 1.08496183
d-24, 1.17440114
d-24, 1.13943939
d-24, &
121 9.71696981
d-25, 7.09593688
d-25, 4.31376399
d-25, 2.12708486
d-25, &
122 8.47429567
d-26, 3.17608104
d-26, 1.95898842
d-26, 1.98064242
d-26, &
123 1.67706555
d-26, 8.99126003
d-27, 3.29773878
d-27, 1.28896127
d-27, &
124 8.51169698
d-28, 7.53520167
d-28, 6.18268143
d-28, 4.30034650
d-28, &
125 2.78152409
d-28, 1.95437088
d-28, 1.65896278
d-28, 1.68740181
d-28, &
126 1.76054383
d-28, 1.63978419
d-28, 1.32880591
d-28, 1.00833205
d-28, &
127 7.82252806
d-29, 6.36181741
d-29, 5.34633869
d-29, 4.58013864
d-29, &
128 3.97833422
d-29, 3.49414760
d-29, 3.09790940
d-29, 2.76786227
d-29, &
129 2.48806269
d-29, 2.24823367
d-29, 2.04016653
d-29, 1.85977413
d-29, &
130 1.70367499
d-29, 1.56966125
d-29, 1.45570643
d-29, 1.35964565
d-29, &
131 1.27879263
d-29, 1.21016980
d-29, 1.15132499
d-29, 1.09959628
d-29, &
132 1.05307482
d-29, 1.01040261
d-29, 9.70657096
d-30, 9.33214234
d-30, &
133 8.97689427
d-30, 8.63761192
d-30, 8.31149879
d-30, 7.85162401
d-30, &
134 7.53828281
d-30, 7.23559452
d-30, 6.94341530
d-30, 6.66137038
d-30, &
135 6.38929156
d-30, 6.12669083
d-30, 5.87346434
d-30, 5.62943622
d-30, &
138 data f_193 / 6.40066486
d-32, 4.92737300
d-31, 2.95342934
d-30, 1.28061594
d-29, &
139 3.47747667
d-29, 5.88554792
d-29, 7.72171179
d-29, 9.75609282
d-29, &
140 1.34318963
d-28, 1.96252638
d-28, 2.70163878
d-28, 3.63192965
d-28, &
141 5.28087341
d-28, 8.37821446
d-28, 1.39089159
d-27, 2.31749718
d-27, &
142 3.77510689
d-27, 5.85198594
d-27, 8.26021568
d-27, 1.04870405
d-26, &
143 1.25209374
d-26, 1.47406787
d-26, 1.77174067
d-26, 2.24098537
d-26, &
144 3.05926105
d-26, 4.50018853
d-26, 6.84720216
d-26, 1.00595861
d-25, &
145 1.30759222
d-25, 1.36481773
d-25, 1.15943558
d-25, 1.01467304
d-25, &
146 1.04092532
d-25, 1.15071251
d-25, 1.27416033
d-25, 1.38463476
d-25, &
147 1.47882726
d-25, 1.57041238
d-25, 1.69786224
d-25, 1.94970397
d-25, &
148 2.50332918
d-25, 3.58321431
d-25, 5.18061550
d-25, 6.60405549
d-25, &
149 6.64085365
d-25, 4.83825816
d-25, 2.40545020
d-25, 8.59534098
d-26, &
150 2.90920638
d-26, 1.33204845
d-26, 9.03933926
d-27, 7.78910836
d-27, &
151 7.29342321
d-27, 7.40267022
d-27, 8.05279981
d-27, 8.13829291
d-27, &
152 6.92634262
d-27, 5.12521880
d-27, 3.59527615
d-27, 2.69617560
d-27, &
153 2.84432713
d-27, 5.06697306
d-27, 1.01281903
d-26, 1.63526978
d-26, &
154 2.06759342
d-26, 2.19482312
d-26, 2.10050611
d-26, 1.89837248
d-26, &
155 1.66347131
d-26, 1.43071097
d-26, 1.21518419
d-26, 1.02078343
d-26, &
156 8.46936184
d-27, 6.93015742
d-27, 5.56973237
d-27, 4.38951754
d-27, &
157 3.38456457
d-27, 2.55309556
d-27, 1.88904224
d-27, 1.38057546
d-27, &
158 1.00718330
d-27, 7.43581116
d-28, 5.63562931
d-28, 4.43359435
d-28, &
159 3.63923535
d-28, 3.11248143
d-28, 2.75586846
d-28, 2.50672237
d-28, &
160 2.32419348
d-28, 2.18325682
d-28, 2.06834486
d-28, 1.93497044
d-28, &
161 1.84540751
d-28, 1.76356504
d-28, 1.68741425
d-28, 1.61566157
d-28, &
162 1.54754523
d-28, 1.48249410
d-28, 1.42020176
d-28, 1.36045230
d-28, &
165 data f_211 / 4.74439912
d-42, 1.95251522
d-40, 5.19700194
d-39, 8.53120166
d-38, &
166 7.72745727
d-37, 4.04158559
d-36, 1.64853511
d-35, 8.56295439
d-35, &
167 1.17529722
d-33, 2.16867729
d-32, 2.40472264
d-31, 1.56418133
d-30, &
168 7.20032889
d-30, 2.65838271
d-29, 8.33196904
d-29, 2.26128236
d-28, &
169 5.24295811
d-28, 9.77791121
d-28, 1.35913489
d-27, 1.43957785
d-27, &
170 1.37591544
d-27, 1.49029886
d-27, 2.06183401
d-27, 3.31440622
d-27, &
171 5.42497318
d-27, 8.41100374
d-27, 1.17941366
d-26, 1.49269794
d-26, &
172 1.71506074
d-26, 1.71266353
d-26, 1.51434781
d-26, 1.36766622
d-26, &
173 1.33483562
d-26, 1.36834518
d-26, 1.45829002
d-26, 1.62575306
d-26, &
174 1.88773347
d-26, 2.22026986
d-26, 2.54930499
d-26, 2.80758138
d-26, &
175 3.06176409
d-26, 3.62799792
d-26, 5.13226109
d-26, 8.46260744
d-26, &
176 1.38486586
d-25, 1.86192535
d-25, 1.78007934
d-25, 1.16548409
d-25, &
177 5.89293257
d-26, 2.69952884
d-26, 1.24891081
d-26, 6.41273176
d-27, &
178 4.08282914
d-27, 3.26463328
d-27, 2.76230280
d-27, 2.08986882
d-27, &
179 1.37658470
d-27, 8.48489381
d-28, 5.19304217
d-28, 3.19312514
d-28, &
180 2.02968197
d-28, 1.50171666
d-28, 1.39164218
d-28, 1.42448821
d-28, &
181 1.41714519
d-28, 1.33341059
d-28, 1.20759270
d-28, 1.07259692
d-28, &
182 9.44895400
d-29, 8.29030041
d-29, 7.25440631
d-29, 6.33479483
d-29, &
183 5.51563757
d-29, 4.79002469
d-29, 4.14990482
d-29, 3.59384972
d-29, &
184 3.12010860
d-29, 2.72624742
d-29, 2.40734791
d-29, 2.15543565
d-29, &
185 1.95921688
d-29, 1.80682882
d-29, 1.68695662
d-29, 1.59020936
d-29, &
186 1.50940886
d-29, 1.43956179
d-29, 1.37731622
d-29, 1.32049043
d-29, &
187 1.26771875
d-29, 1.21803879
d-29, 1.17074716
d-29, 1.10507836
d-29, &
188 1.06022834
d-29, 1.01703080
d-29, 9.75436986
d-30, 9.35349257
d-30, &
189 8.96744546
d-30, 8.59527489
d-30, 8.23678940
d-30, 7.89144480
d-30, &
192 data f_304 / 3.62695850
d-32, 2.79969087
d-31, 1.68340584
d-30, 7.32681440
d-30, &
193 1.99967770
d-29, 3.41296785
d-29, 4.55409104
d-29, 5.94994635
d-29, &
194 8.59864963
d-29, 1.39787633
d-28, 3.17701965
d-28, 1.14474920
d-27, &
195 4.44845958
d-27, 1.54785841
d-26, 4.70265345
d-26, 1.24524365
d-25, &
196 2.81535352
d-25, 5.10093666
d-25, 6.83545307
d-25, 6.82110329
d-25, &
197 5.66886188
d-25, 4.36205513
d-25, 3.29265688
d-25, 2.49802368
d-25, &
198 1.92527113
d-25, 1.51058572
d-25, 1.20596047
d-25, 9.76884267
d-26, &
199 7.89979266
d-26, 6.18224289
d-26, 4.67298332
d-26, 3.57934505
d-26, &
200 2.84535785
d-26, 2.32853022
d-26, 1.95228514
d-26, 1.67880071
d-26, &
201 1.47608785
d-26, 1.32199691
d-26, 1.20070960
d-26, 1.09378177
d-26, &
202 1.00031730
d-26, 9.62434001
d-27, 1.05063954
d-26, 1.27267143
d-26, &
203 1.45923057
d-26, 1.36746707
d-26, 1.03466970
d-26, 6.97647829
d-27, &
204 4.63141039
d-27, 3.19031994
d-27, 2.33373613
d-27, 1.81589079
d-27, &
205 1.48446917
d-27, 1.26611478
d-27, 1.12617468
d-27, 1.03625148
d-27, &
206 9.61400595
d-28, 8.79016231
d-28, 7.82612130
d-28, 6.73762960
d-28, &
207 5.59717956
d-28, 4.53010243
d-28, 3.65712196
d-28, 3.00958686
d-28, &
208 2.54011502
d-28, 2.18102277
d-28, 1.88736437
d-28, 1.63817539
d-28, &
209 1.42283147
d-28, 1.23631916
d-28, 1.07526003
d-28, 9.36797928
d-29, &
210 8.18565660
d-29, 7.18152734
d-29, 6.32523238
d-29, 5.59513985
d-29, &
211 4.96614048
d-29, 4.42518826
d-29, 3.95487628
d-29, 3.54690294
d-29, &
212 3.18953930
d-29, 2.87720933
d-29, 2.60186750
d-29, 2.36011522
d-29, &
213 2.14717806
d-29, 1.95905217
d-29, 1.79287981
d-29, 1.64562262
d-29, &
214 1.51489425
d-29, 1.39876064
d-29, 1.29496850
d-29, 1.18665438
d-29, &
215 1.10240474
d-29, 1.02643099
d-29, 9.57780996
d-30, 8.95465151
d-30, &
216 8.38950190
d-30, 7.87283711
d-30, 7.40136507
d-30, 6.96804279
d-30, &
219 data f_335 / 2.46882661
d-32, 1.89476632
d-31, 1.13216502
d-30, 4.89532008
d-30, &
220 1.32745970
d-29, 2.25390335
d-29, 3.00511672
d-29, 3.96035934
d-29, &
221 5.77977656
d-29, 8.58600736
d-29, 1.14083000
d-28, 1.48644411
d-28, &
222 2.15788823
d-28, 3.51628877
d-28, 6.12200698
d-28, 1.08184987
d-27, &
223 1.85590697
d-27, 2.91679107
d-27, 3.94405396
d-27, 4.63610680
d-27, &
224 5.13824456
d-27, 5.66602209
d-27, 6.30009232
d-27, 7.03422868
d-27, &
225 7.77973918
d-27, 8.32371831
d-27, 8.56724316
d-27, 8.62601374
d-27, &
226 8.13308844
d-27, 6.53188216
d-27, 4.55197029
d-27, 3.57590087
d-27, &
227 3.59571707
d-27, 4.03502770
d-27, 4.54366411
d-27, 4.96914990
d-27, &
228 5.24601170
d-27, 5.39979250
d-27, 5.43023669
d-27, 5.26235042
d-27, &
229 4.91585495
d-27, 4.52628362
d-27, 4.13385020
d-27, 3.67538967
d-27, &
230 3.39939742
d-27, 3.81284533
d-27, 5.02332701
d-27, 6.19438602
d-27, &
231 6.49613071
d-27, 6.04010475
d-27, 5.24664275
d-27, 4.37225997
d-27, &
232 3.52957182
d-27, 2.76212276
d-27, 2.08473158
d-27, 1.50850518
d-27, &
233 1.04602472
d-27, 7.13091243
d-28, 5.34289645
d-28, 5.21079581
d-28, &
234 6.22246365
d-28, 6.99555864
d-28, 6.29665489
d-28, 4.45077026
d-28, &
235 2.67046793
d-28, 1.52774686
d-28, 9.18061770
d-29, 6.09116074
d-29, &
236 4.48562572
d-29, 3.59463696
d-29, 3.05820218
d-29, 2.70766652
d-29, &
237 2.46144034
d-29, 2.27758450
d-29, 2.13331183
d-29, 2.01537836
d-29, &
238 1.91566180
d-29, 1.82893912
d-29, 1.75167748
d-29, 1.68136168
d-29, &
239 1.61615595
d-29, 1.55481846
d-29, 1.49643236
d-29, 1.44046656
d-29, &
240 1.38657085
d-29, 1.33459068
d-29, 1.28447380
d-29, 1.23615682
d-29, &
241 1.18963296
d-29, 1.14478976
d-29, 1.10146637
d-29, 1.04039479
d-29, &
242 9.98611410
d-30, 9.58205147
d-30, 9.19202009
d-30, 8.81551313
d-30, &
243 8.45252127
d-30, 8.10224764
d-30, 7.76469090
d-30, 7.43954323
d-30, &
249 data t_iris / 4. , 4.1 , 4.2 , 4.3 , 4.40000001, &
250 4.50000001, 4.60000001, 4.70000001, 4.80000001, 4.90000001, &
251 5.00000001, 5.10000002, 5.20000002, 5.30000002, 5.40000002, &
252 5.50000002, 5.60000002, 5.70000003, 5.80000003, 5.90000003, &
253 6.00000003, 6.10000003, 6.20000003, 6.30000003, 6.40000004, &
254 6.50000004, 6.60000004, 6.70000004, 6.80000004, 6.90000004, &
255 7.00000004, 7.10000005, 7.20000005, 7.30000005, 7.40000005, &
256 7.50000005, 7.60000005, 7.70000006, 7.80000006, 7.90000006, &
259 data f_1354 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
260 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
261 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
262 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
263 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
264 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 1.09503647
d-39, &
265 5.47214550
d-36, 2.42433983
d-33, 2.75295034
d-31, 1.21929718
d-29, &
266 2.48392125
d-28, 2.33268145
d-27, 8.68623633
d-27, 1.00166284
d-26, &
267 3.63126633
d-27, 7.45174807
d-28, 1.38224064
d-28, 2.69270994
d-29, &
268 5.53314977
d-30, 1.15313092
d-30, 2.34195788
d-31, 4.48242942
d-32, &
274 data t_eis1 / 1.99526231
d+05, 2.23872114
d+05, 2.51188643
d+05, 2.81838293
d+05, &
275 3.16227766
d+05, 3.54813389
d+05, 3.98107171
d+05, 4.46683592
d+05, &
276 5.01187234
d+05, 5.62341325
d+05, 6.30957344
d+05, 7.07945784
d+05, &
277 7.94328235
d+05, 8.91250938
d+05, 1.00000000
d+06, 1.12201845
d+06, &
278 1.25892541
d+06, 1.41253754
d+06, 1.58489319
d+06, 1.77827941
d+06, &
279 1.99526231
d+06, 2.23872114
d+06, 2.51188643
d+06, 2.81838293
d+06, &
280 3.16227766
d+06, 3.54813389
d+06, 3.98107171
d+06, 4.46683592
d+06, &
281 5.01187234
d+06, 5.62341325
d+06, 6.30957344
d+06, 7.07945784
d+06, &
282 7.94328235
d+06, 8.91250938
d+06, 1.00000000
d+07, 1.12201845
d+07, &
283 1.25892541
d+07, 1.41253754
d+07, 1.58489319
d+07, 1.77827941
d+07, &
284 1.99526231
d+07, 2.23872114
d+07, 2.51188643
d+07, 2.81838293
d+07, &
285 3.16227766
d+07, 3.54813389
d+07, 3.98107171
d+07, 4.46683592
d+07, &
286 5.01187234
d+07, 5.62341325
d+07, 6.30957344
d+07, 7.07945784
d+07, &
287 7.94328235
d+07, 8.91250938
d+07, 1.00000000
d+08, 1.12201845
d+08, &
288 1.25892541
d+08, 1.41253754
d+08, 1.58489319
d+08, 1.77827941
d+08 /
290 data t_eis2 / 1.99526231
d+06, 2.23872114
d+06, 2.51188643
d+06, 2.81838293
d+06, &
291 3.16227766
d+06, 3.54813389
d+06, 3.98107171
d+06, 4.46683592
d+06, &
292 5.01187234
d+06, 5.62341325
d+06, 6.30957344
d+06, 7.07945784
d+06, &
293 7.94328235
d+06, 8.91250938
d+06, 1.00000000
d+07, 1.12201845
d+07, &
294 1.25892541
d+07, 1.41253754
d+07, 1.58489319
d+07, 1.77827941
d+07, &
295 1.99526231
d+07, 2.23872114
d+07, 2.51188643
d+07, 2.81838293
d+07, &
296 3.16227766
d+07, 3.54813389
d+07, 3.98107171
d+07, 4.46683592
d+07, &
297 5.01187234
d+07, 5.62341325
d+07, 6.30957344
d+07, 7.07945784
d+07, &
298 7.94328235
d+07, 8.91250938
d+07, 1.00000000
d+08, 1.12201845
d+08, &
299 1.25892541
d+08, 1.41253754
d+08, 1.58489319
d+08, 1.77827941
d+08, &
300 1.99526231
d+08, 2.23872114
d+08, 2.51188643
d+08, 2.81838293
d+08, &
301 3.16227766
d+08, 3.54813389
d+08, 3.98107171
d+08, 4.46683592
d+08, &
302 5.01187234
d+08, 5.62341325
d+08, 6.30957344
d+08, 7.07945784
d+08, &
303 7.94328235
d+08, 8.91250938
d+08, 1.00000000
d+09, 1.12201845
d+09, &
304 1.25892541
d+09, 1.41253754
d+09, 1.58489319
d+09, 1.77827941
d+09 /
306 data f_263 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
307 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
308 0.00000000
d+00, 4.46454917
d-45, 3.26774829
d-42, 1.25292566
d-39, &
309 2.66922338
d-37, 3.28497742
d-35, 2.38677554
d-33, 1.03937729
d-31, &
310 2.75075687
d-30, 4.47961733
d-29, 4.46729177
d-28, 2.64862689
d-27, &
311 8.90863800
d-27, 1.72437548
d-26, 2.22217752
d-26, 2.27999477
d-26, &
312 2.08264363
d-26, 1.78226687
d-26, 1.45821699
d-26, 1.14675379
d-26, &
313 8.63082492
d-27, 6.15925429
d-27, 4.11252514
d-27, 2.51530564
d-27, &
314 1.37090986
d-27, 6.42443134
d-28, 2.48392636
d-28, 7.59187874
d-29, &
315 1.77852938
d-29, 3.23945221
d-30, 4.90533903
d-31, 6.75458158
d-32, &
316 9.06878868
d-33, 1.23927474
d-33, 1.75769395
d-34, 2.60710914
d-35, &
317 4.04318030
d-36, 6.53500581
d-37, 1.09365022
d-37, 1.88383322
d-38, &
318 3.31425233
d-39, 5.90964084
d-40, 1.06147549
d-40, 1.90706170
d-41, &
319 3.41331584
d-42, 6.07310718
d-43, 1.07364738
d-43, 1.89085498
d-44, &
320 3.32598922
d-45, 5.87125640
d-46, 0.00000000
d+00, 0.00000000
d+00 /
322 data f_264 / 0.00000000
d+00, 2.81670057
d-46, 1.28007268
d-43, 2.54586603
d-41, &
323 2.67887256
d-39, 1.68413285
d-37, 6.85702304
d-36, 1.91797284
d-34, &
324 3.84675839
d-33, 5.69939170
d-32, 6.36224608
d-31, 5.39176489
d-30, &
325 3.45478458
d-29, 1.64848693
d-28, 5.71476364
d-28, 1.39909997
d-27, &
326 2.37743056
d-27, 2.86712530
d-27, 2.65206348
d-27, 2.07175767
d-27, &
327 1.47866767
d-27, 1.01087374
d-27, 6.79605811
d-28, 4.54746770
d-28, &
328 3.04351751
d-28, 2.03639149
d-28, 1.35940991
d-28, 9.01451939
d-29, &
329 5.91289972
d-29, 3.81821178
d-29, 2.41434696
d-29, 1.48871220
d-29, &
330 8.93362094
d-30, 5.21097445
d-30, 2.95964719
d-30, 1.64278748
d-30, &
331 8.95571660
d-31, 4.82096011
d-31, 2.57390991
d-31, 1.36821781
d-31, &
332 7.27136350
d-32, 3.87019426
d-32, 2.06883430
d-32, 1.11228884
d-32, &
333 6.01883313
d-33, 3.27790676
d-33, 1.79805012
d-33, 9.93085346
d-34, &
334 5.52139556
d-34, 3.08881387
d-34, 1.73890315
d-34, 9.84434964
d-35, &
335 5.60603378
d-35, 3.20626492
d-35, 1.84111068
d-35, 0.00000000
d+00, &
336 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00 /
338 data f_192 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 4.35772105
d-44, &
339 1.26162319
d-41, 1.97471205
d-39, 1.83409019
d-37, 1.08206288
d-35, &
340 4.27914363
d-34, 1.17943846
d-32, 2.32565755
d-31, 3.33087991
d-30, &
341 3.47013260
d-29, 2.60375866
d-28, 1.37737127
d-27, 5.01053913
d-27, &
342 1.23479810
d-26, 2.11310542
d-26, 2.71831513
d-26, 2.89851163
d-26, &
343 2.77312376
d-26, 2.50025229
d-26, 2.18323661
d-26, 1.86980322
d-26, &
344 1.58035034
d-26, 1.31985651
d-26, 1.08733133
d-26, 8.81804906
d-27, &
345 7.00417973
d-27, 5.43356567
d-27, 4.09857884
d-27, 2.99651764
d-27, &
346 2.11902962
d-27, 1.45014127
d-27, 9.62291023
d-28, 6.21548647
d-28, &
347 3.92807578
d-28, 2.44230375
d-28, 1.50167782
d-28, 9.17611405
d-29, &
348 5.58707641
d-29, 3.40570915
d-29, 2.08030862
d-29, 1.27588676
d-29, &
349 7.86535588
d-30, 4.87646151
d-30, 3.03888897
d-30, 1.90578649
d-30, &
350 1.20195947
d-30, 7.61955060
d-31, 4.85602199
d-31, 3.11049969
d-31, &
351 2.00087065
d-31, 1.29223740
d-31, 8.37422008
d-32, 0.00000000
d+00, &
352 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00 /
354 data f_255 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 1.76014287
d-44, &
355 5.07057938
d-42, 7.90473970
d-40, 7.31852999
d-38, 4.30709255
d-36, &
356 1.70009061
d-34, 4.67925160
d-33, 9.21703546
d-32, 1.31918676
d-30, &
357 1.37393161
d-29, 1.03102379
d-28, 5.45694018
d-28, 1.98699648
d-27, &
358 4.90346776
d-27, 8.40524725
d-27, 1.08321456
d-26, 1.15714525
d-26, &
359 1.10905152
d-26, 1.00155023
d-26, 8.75799161
d-27, 7.50935839
d-27, &
360 6.35253533
d-27, 5.30919268
d-27, 4.37669455
d-27, 3.55185164
d-27, &
361 2.82347055
d-27, 2.19257595
d-27, 1.65589541
d-27, 1.21224987
d-27, &
362 8.58395132
d-28, 5.88163935
d-28, 3.90721447
d-28, 2.52593407
d-28, &
363 1.59739995
d-28, 9.93802874
d-29, 6.11343388
d-29, 3.73711135
d-29, &
364 2.27618743
d-29, 1.38793199
d-29, 8.48060787
d-30, 5.20305940
d-30, &
365 3.20867365
d-30, 1.99011277
d-30, 1.24064551
d-30, 7.78310544
d-31, &
366 4.91013681
d-31, 3.11338381
d-31, 1.98451675
d-31, 1.27135460
d-31, &
367 8.17917486
d-32, 5.28280497
d-32, 3.42357159
d-32, 0.00000000
d+00, &
368 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00 /
373 integer,
intent(in) :: ixI^L, ixO^L
374 double precision,
intent(in) :: w(ixI^S,nw)
375 double precision,
intent(in) :: x(ixI^S,1:ndim)
376 double precision,
intent(out):: res(ixI^S)
384 integer,
intent(in) :: ixI^L, ixO^L
385 double precision,
intent(in) :: w(ixI^S, nw)
386 double precision,
intent(out):: val1(ixI^S), val2(ixI^S)
392 procedure(
get_subr1),
pointer,
nopass :: get_rho => null()
393 procedure(
get_subr1),
pointer,
nopass :: get_pthermal => null()
394 procedure(
get_subr1),
pointer,
nopass :: get_var_rfactor => null()
405 double precision,
allocatable :: source(:^d&)
406 double precision,
allocatable :: opacity(:^d&)
407 double precision,
allocatable :: sourcev(:^d&)
408 double precision,
allocatable :: xface1(:),xface2(:),xface3(:)
409 double precision,
allocatable :: rface(:),thetaface(:),phiface(:)
410 double precision,
allocatable :: rface2(:),theta_cos(:),phi_sin(:),phi_cos(:)
411 double precision :: box_min(1:3)=0.d0
412 double precision :: box_max(1:3)=0.d0
417 logical :: has_pixels=.false.
429 character(len=*),
intent(in) :: datatype
432 call mpistop(
"bad radiation_transfer")
448 case(
'cart',
'cart_dda')
450 case(
'spherical',
'sph_intersection')
452 case(
'sph_dda',
'spherical_dda')
462 call mpistop(
"bad emission_model")
466 call mpistop(
"tau and absorption-fraction output need thick transfer")
470 call mpistop(
"radsyn_pixel_batch must be positive")
473 call mpistop(
"radsyn_segment_batch_factor must be non-negative")
476 call mpistop(
"radsyn_segment_memory_mb must be positive")
479 call mpistop(
"radsyn_segment_comm_factor must be positive")
484 call mpistop(
"instrument_postprocess currently needs dat-resolution EUV images")
487 call mpistop(
"instrument_postprocess is not yet supported for spherical rays")
490 call mpistop(
"instrument_postprocess currently supports only EUV AIA or radio_ff images")
493 call mpistop(
"radio_ff instrument_postprocess needs radio_beam_fwhm > 0 arcsec")
501 if (datatype /=
'image_euv' .and. datatype /=
'spectrum_euv')
then
502 call mpistop(
"emission_model=euv_aia is only valid for EUV synthesis")
505 if (datatype /=
'image_whitelight')
then
506 call mpistop(
"emission_model=white_light is only valid for white-light synthesis")
509 if (datatype /=
'image_euv')
then
510 call mpistop(
"emission_model=radio_ff currently reuses EUV-image convert types")
513 call mpistop(
"emission_model=radio_ff needs radio_frequency > 0")
515 case(
'pseudo_current')
516 if (datatype /=
'image_euv')
then
517 call mpistop(
"emission_model=pseudo_current is only valid for EUV-image convert types")
520 call mpistop(
"emission_model=pseudo_current currently supports only thin transfer")
526 call mpistop(
"ray_method=cart needs Cartesian EUV slab images")
531 call mpistop(
"ray_method=spherical currently needs 3D spherical grids")
534 call mpistop(
"ray_method=spherical currently needs 3D spherical grids")
539 call mpistop(
"bad ray_method=spherical mode")
542 call mpistop(
"ray_method=spherical currently supports only EUV AIA emission")
544 if (xprobmin2<=1.
d-10 .or. xprobmax2>=dpi-1.
d-10)
then
545 call mpistop(
"ray_method=spherical does not support polar-axis crossing domains")
547 if (xprobmax3<=xprobmin3 .or. xprobmax3-xprobmin3>=2.d0*dpi-1.
d-10)
then
548 call mpistop(
"ray_method=spherical does not support phi-wrapping domains")
553 if (datatype /=
'image_euv')
then
554 call mpistop(
"thick transfer is only defined for EUV images")
559 if (.not.
slab)
call mpistop(
"cartesian thick EUV currently needs slab output")
561 call mpistop(
"thick EUV currently needs Cartesian dat_resolution output")
567 call mpistop(
"thick EUV currently needs x/y/z-aligned LOS")
574 double precision,
intent(in) :: emissivity,opacity,path_length
575 double precision,
intent(inout) :: intensity,tau
577 double precision :: dtau
579 if (path_length<=zero)
return
581 dtau=max(zero,opacity)*path_length
587 trim(emission_model)/=
'radio_ff' .and. &
592 logical,
intent(in) :: has_doppler,has_thick
595 if (has_doppler) num_outputs=num_outputs+1
596 if (has_thick .and. output_tau) num_outputs=num_outputs+1
597 if (has_thick .and. output_absorption_fraction) num_outputs=num_outputs+1
601 integer,
intent(in) :: nI1,nI2
602 double precision,
intent(in) :: EUV(nI1,nI2),unitv
603 double precision,
intent(inout) :: Dpl(nI1,nI2)
609 if (euv(ix1,ix2)/=zero)
then
610 dpl(ix1,ix2)=(dpl(ix1,ix2)/euv(ix1,ix2))*unitv
614 if (abs(dpl(ix1,ix2))<smalldouble) dpl(ix1,ix2)=zero
620 integer,
intent(in) :: nI1,nI2
621 double precision,
intent(in) :: EUV(nI1,nI2),EUVthin(nI1,nI2),smallflux
622 double precision,
intent(out) :: Absorption(nI1,nI2)
623 logical,
intent(in),
optional :: cap_to_one
626 logical :: cap_absorption
629 cap_absorption=.false.
630 if (
present(cap_to_one)) cap_absorption=cap_to_one
633 if (euvthin(ix1,ix2)>smallflux)
then
634 absorption(ix1,ix2)=max(zero,(euvthin(ix1,ix2)-euv(ix1,ix2))/euvthin(ix1,ix2))
635 if (cap_absorption) absorption(ix1,ix2)=min(one,absorption(ix1,ix2))
641 subroutine pack_euv_image_outputs(nI1,nI2,EUV,wI,smallflux,has_doppler,has_thick,Dpl,Tau,EUVthin,&
643 integer,
intent(in) :: nI1,nI2
644 double precision,
intent(in) :: EUV(nI1,nI2),smallflux
645 double precision,
intent(inout) :: wI(:,:,:)
646 logical,
intent(in) :: has_doppler,has_thick
647 double precision,
intent(in),
optional :: Dpl(nI1,nI2),Tau(nI1,nI2),EUVthin(nI1,nI2)
648 logical,
intent(in),
optional :: cap_absorption
651 double precision,
allocatable :: Absorption(:,:)
656 if (has_doppler)
then
657 if (.not.
present(dpl))
call mpistop(
"Doppler output requested without Doppler image")
661 if (has_thick .and. output_tau)
then
662 if (.not.
present(tau))
call mpistop(
"tau output requested without tau image")
666 if (has_thick .and. output_absorption_fraction)
then
667 if (.not.
present(euvthin))
call mpistop(
"absorption output requested without thin image")
668 allocate(absorption(ni1,ni2))
671 wi(:,:,iw)=absorption(:,:)
672 deallocate(absorption)
677 integer,
intent(out) :: pixel_batch_target,segment_batch_target,segment_comm_target
679 pixel_batch_target=max(1,radsyn_pixel_batch)
680 if (radsyn_segment_batch_factor>0)
then
681 segment_batch_target=max(128,radsyn_segment_batch_factor*pixel_batch_target)
683 segment_batch_target=max(128,int(min(dble(huge(segment_batch_target)),&
684 max(128.d0,radsyn_segment_memory_mb*1048576.d0/256.d0))))
686 segment_comm_target=max(128,radsyn_segment_comm_factor*pixel_batch_target)
690 double precision,
intent(in) :: tau
700 double precision,
intent(in) :: argument
702 if (argument<-700.d0)
then
704 else if (argument>700.d0)
then
712 double precision,
intent(in) :: exponent
714 if (exponent>300.d0)
then
716 else if (exponent<-300.d0)
then
724 double precision,
intent(in) :: temperature
725 integer,
intent(in) :: n_table
726 double precision,
intent(in) :: t_table(n_table),f_table(n_table)
727 logical,
intent(in) :: log_temperature,log_response
729 integer :: ilo,ihi,imid
730 double precision :: temp_lookup,response_lookup,flo,fhi
733 if (temperature<=zero)
return
734 if (log_temperature)
then
735 temp_lookup=log10(temperature)
737 temp_lookup=temperature
739 if (temp_lookup<t_table(1) .or. temp_lookup>t_table(n_table))
return
740 if (temp_lookup==t_table(n_table))
then
741 if (log_response)
then
742 response_lookup=log10(max(f_table(n_table),1.d-99))
744 response_lookup=f_table(n_table)
751 if (temp_lookup>=t_table(imid))
then
757 if (log_response)
then
758 flo=log10(max(f_table(ilo),1.d-99))
759 fhi=log10(max(f_table(ilo+1),1.d-99))
764 response_lookup=flo*(temp_lookup-t_table(ilo+1))/(t_table(ilo)-t_table(ilo+1))+&
765 fhi*(temp_lookup-t_table(ilo))/(t_table(ilo+1)-t_table(ilo))
768 if (log_response)
then
777 integer,
intent(in) :: ixI^L, ixO^L, n_table
778 double precision,
intent(in) :: Te(ixI^S),t_table(n_table),f_table(n_table)
779 double precision,
intent(inout) :: flux(ixI^S)
780 logical,
intent(in) :: log_temperature,log_response
783 double precision :: GT
785 {
do ix^db=ixomin^db,ixomax^db\}
787 flux(ix^d)=flux(ix^d)*gt
788 if (flux(ix^d)<zero) flux(ix^d)=zero
797 integer,
intent(in) :: wl
798 integer,
intent(in) :: ixI^L, ixO^L
799 double precision,
intent(in) :: x(ixI^S,1:ndim)
800 double precision,
intent(in) :: w(ixI^S,1:nw)
802 double precision,
intent(out) :: kappa(ixI^S)
805 double precision :: pth(ixI^S),Te(ixI^S),Ne(ixI^S)
806 double precision :: wave_ratio,s_H1,s_He1,s_He2,Pe
807 double precision :: log_H21,log_He21,log_He32,log_He321,logScaleHe,w_H21
808 double precision :: term0,term1,term2,denHe,i0,j1,j2,be
809 double precision :: N_H,N_H1,N_He1,N_He2
810 double precision,
parameter :: Xe_H21=13.6d0, xe_he21=24.587d0, xe_he32=54.416d0
811 double precision,
parameter :: rHe=0.1d0
812 double precision,
parameter :: sigma_H1=5.16d-20, sigma_he1=9.25d-19, sigma_he2=7.17d-19
814 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
815 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
816 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
817 te(ixo^s)=pth(ixo^s)/(ne(ixo^s)*te(ixo^s))*unit_temperature
819 double precision :: nH_dummy(ixI^S)
820 call eos%get_ne_nH(ixi^l, ixo^l, w, ne, nh_dummy)
823 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
825 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
828 wave_ratio=dble(wl)/171.d0
832 if (wl<=912) s_h1=wave_ratio**3*sigma_h1
833 if (wl<=504) s_he1=wave_ratio**2*sigma_he1
834 if (wl<=228) s_he2=wave_ratio**2.75d0*sigma_he2
837 {
do ix^db=ixomin^db,ixomax^db\}
838 if (te(ix^d)>zero .and. ne(ix^d)>zero)
then
839 pe=ne(ix^d)*kb_cgs*te(ix^d)
841 log_h21=2.5d0*log10(te(ix^d))-5040.d0*xe_h21/te(ix^d)-log10(pe)-0.48d0
842 log_he21=log10(4.d0)+2.5d0*log10(te(ix^d))-5040.d0*xe_he21/te(ix^d)-log10(pe)-0.48d0
843 log_he32=2.5d0*log10(te(ix^d))-5040.d0*xe_he32/te(ix^d)-log10(pe)-0.48d0
845 i0=w_h21/(1.d0+w_h21)
846 log_he321=log_he21+log_he32
847 logscalehe=max(zero,log_he21,log_he321)
851 denhe=term0+term1+term2
859 be=i0+rhe*(j1+2.d0*j2)
860 if (be>smalldouble)
then
863 n_he1=(1.d0-j1-j2)*rhe*n_h
865 kappa(ix^d)=max(zero,n_h1*s_h1+n_he1*s_he1+n_he2*s_he2)
873 integer,
intent(in) :: igrid
874 integer,
intent(in) :: ixI^L, ixO^L
875 double precision,
intent(in) :: w(ixI^S,1:nw)
876 double precision,
intent(out) :: source(ixI^S)
878 integer :: ix^D,idir,idirmin,idirmin0
879 double precision :: current(ixI^S,7-2*ndir:3)
881 if (.not.
allocated(iw_mag))
then
882 call mpistop(
"emission_model=pseudo_current needs magnetic-field variables")
887 call curlvector(w(ixi^s,iw_mag(1:ndir)),ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
889 current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+ps(igrid)%J0(ixo^s,idirmin0:3)
893 {
do ix^db=ixomin^db,ixomax^db\}
895 source(ix^d)=source(ix^d)+current(ix^d,idir)**2
903 integer,
intent(in) :: ixI^L, ixO^L
904 double precision,
intent(in) :: x(ixI^S,1:ndim)
905 double precision,
intent(in) :: w(ixI^S,1:nw)
907 double precision,
intent(out) :: source(ixI^S),kappa(ixI^S)
910 double precision :: pth(ixI^S),Te(ixI^S),Ne(ixI^S)
911 double precision :: nH_dummy(ixI^S),gff
913 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
914 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
915 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
916 te(ixo^s)=pth(ixo^s)/(ne(ixo^s)*te(ixo^s))*unit_temperature
917 call eos%get_ne_nH(ixi^l,ixo^l,w,ne,nh_dummy)
919 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
921 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
926 {
do ix^db=ixomin^db,ixomax^db\}
927 if (te(ix^d)>zero .and. ne(ix^d)>zero)
then
928 if (te(ix^d)<2.d5)
then
929 gff=18.2d0+1.5d0*log(te(ix^d))-log(radio_frequency)
931 gff=24.5d0+log(te(ix^d))-log(radio_frequency)
934 kappa(ix^d)=9.78d-3*ne(ix^d)**2*gff/(radio_frequency**2*te(ix^d)**1.5d0)
935 source(ix^d)=te(ix^d)*kappa(ix^d)
940 subroutine get_line_info(wl,ion,mass,logTe,line_center,spatial_px,spectral_px,sigma_PSF,width_slit)
952 integer,
intent(in) :: wl
953 integer,
intent(out) :: mass
954 character(len=30),
intent(out) :: ion
955 double precision,
intent(out) :: logTe,line_center,spatial_px,spectral_px
956 double precision,
intent(out) :: sigma_PSF,width_slit
1026 line_center=1354.1d0
1028 spectral_px=12.98
d-3
1035 line_center=262.976d0
1044 line_center=263.765d0
1053 line_center=192.028d0
1062 line_center=255.113d0
1068 call mpistop(
"No information about this line")
1082 integer,
intent(in) :: wl
1083 integer,
intent(in) :: ixI^L, ixO^L
1084 double precision,
intent(in) :: x(ixI^S,1:ndim)
1085 double precision,
intent(in) :: w(ixI^S,1:nw)
1087 double precision,
intent(out) :: flux(ixI^S)
1090 double precision :: pth(ixI^S),Te(ixI^S),Ne(ixI^S)
1092 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
1093 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
1094 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
1098 double precision :: nH_dummy(ixI^S)
1099 call eos%get_ne_nH(ixi^l, ixo^l, w, ne, nh_dummy)
1102 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
1103 flux(ixo^s)=ne(ixo^s)**2
1105 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
1106 flux(ixo^s)=ne(ixo^s)**2
1135 call mpistop(
"Unknown wavelength")
1147 integer,
intent(in) :: ixI^L,ixO^L
1148 integer,
intent(in) :: El,Eu
1149 double precision,
intent(in) :: x(ixI^S,1:ndim)
1150 double precision,
intent(in) :: w(ixI^S,nw)
1152 double precision,
intent(out) :: flux(ixI^S)
1154 integer :: ix^D,ixO^D
1156 double precision :: I0,kb,keV,dE,Ei
1157 double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
1158 double precision :: Ne(ixI^S),gff(ixI^S),fi(ixI^S)
1159 double precision :: EM(ixI^S)
1165 nume=floor((eu-el)/de)
1166 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
1167 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
1168 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
1172 double precision :: nH_dummy(ixI^S)
1173 call eos%get_ne_nH(ixi^l, ixo^l, w, ne, nh_dummy)
1176 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
1177 em(ixo^s)=(ne(ixo^s))**2*1.d6
1179 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
1180 em(ixo^s)=(ne(ixo^s))**2
1182 kbt(ixo^s)=kb*te(ixo^s)/kev
1187 {
do ix^db=ixomin^db,ixomax^db\}
1188 if (kbt(ix^d)>0.01*ei)
then
1189 if(kbt(ix^d)<ei) gff(ix^d)=(kbt(ix^d)/ei)**0.4
1190 fi(ix^d)=(em(ix^d)*gff(ix^d))*
exp_clamped(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
1195 flux(ixo^s)=flux(ixo^s)+fi(ixo^s)*de
1197 flux(ixo^s)=flux(ixo^s)*i0
1204 double precision,
intent(in) :: xbox^L
1206 double precision,
intent(out) :: eflux
1208 double precision :: dxb^D,xb^L
1209 integer :: iigrid,igrid,j
1210 integer :: ixO^L,ixI^L,ix^D
1211 double precision :: eflux_grid,eflux_pe
1213 ^d&iximin^d=
ixglo^d;
1214 ^d&iximax^d=
ixghi^d;
1215 ^d&ixomin^d=ixmlo^d;
1216 ^d&ixomax^d=ixmhi^d;
1218 do iigrid=1,igridstail; igrid=igrids(iigrid);
1222 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)
1223 eflux_pe=eflux_pe+eflux_grid
1225 call mpi_allreduce(eflux_pe,eflux,1,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
1232 integer,
intent(in) :: ixI^L,ixO^L
1233 double precision,
intent(in) :: x(ixI^S,1:ndim),dV(ixI^S)
1234 double precision,
intent(in) :: w(ixI^S,nw)
1235 double precision,
intent(in) :: xbox^L,xb^L
1237 double precision,
intent(out) :: eflux_grid
1239 integer :: ix^D,ixO^D,ixb^L
1240 integer :: iE,numE,j,inbox
1241 double precision :: I0,kb,keV,dE,Ei,El,Eu,A_cgs
1242 double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
1243 double precision :: Ne(ixI^S),EM(ixI^S)
1244 double precision :: gff,fi,erg_SI
1248 {
if (xbmin^d<xboxmax^d .and. xbmax^d>xboxmin^d) inbox=inbox+1\}
1250 if (inbox==ndim)
then
1252 ^d&ixbmin^d=ixomin^d;
1253 ^d&ixbmax^d=ixomax^d;
1254 {
if (xbmax^d>xboxmax^d) ixbmax^d=ixomax^d-ceiling((xbmax^d-xboxmax^d)/
dxlevel(^d))\}
1255 {
if (xbmin^d<xboxmin^d) ixbmin^d=ceiling((xboxmin^d-xbmin^d)/
dxlevel(^d))+ixomin^d\}
1262 el=const_h*const_c/(8.d0*a_cgs)/kev
1263 eu=const_h*const_c/(1.d0*a_cgs)/kev
1265 nume=floor((eu-el)/de)
1266 call fl%get_pthermal(w,x,ixi^l,ixb^l,pth)
1267 call fl%get_rho(w,x,ixi^l,ixb^l,ne)
1268 call fl%get_var_Rfactor(w,x,ixi^l,ixb^l,te)
1272 double precision :: nH_dummy(ixI^S)
1273 call eos%get_ne_nH(ixi^l, ixb^l, w, ne, nh_dummy)
1276 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
1277 em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*(unit_length*1.d2)**3
1279 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
1280 em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*unit_length**3
1282 kbt(ixb^s)=kb*te(ixb^s)/kev
1287 {
do ix^db=ixbmin^db,ixbmax^db\}
1288 if (kbt(ix^d)>1.d-2*ei)
then
1289 if(kbt(ix^d)<ei)
then
1290 gff=(kbt(ix^d)/ei)**0.4
1294 fi=(em(ix^d)*gff)*
exp_clamped(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
1295 eflux_grid=eflux_grid+fi*de*ei
1299 eflux_grid=eflux_grid*kev*erg_si
1308 integer,
intent(in) :: qunit
1310 character(20) :: datatype
1313 character (30) :: ion
1314 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1315 double precision :: xslit,arcsec
1317 datatype=
'spectrum_euv'
1322 if (
mype==0) print *,
'###################################################'
1325 if (
mype==0) print *,
'Systhesizing EUV spectrum (observed by IRIS).'
1326 case (263,264,192,255)
1327 if (
mype==0) print *,
'Systhesizing EUV spectrum (observed by Hinode/EIS).'
1329 call mpistop(
'Wrong wavelength!')
1333 call mpistop(
'Wrong spectrum window!')
1336 if (
mype==0)
write(*,
'(a,f8.3,a)')
' Wavelength: ',linecent,
' Angstrom'
1337 if (
mype==0) print *,
'Unit of EUV flux: DN s^-1 pixel^-1'
1341 write(*,
'(a,f5.3,a,f5.1,a)')
' Supposed pixel: ',wlrsl,
' Angstrom x ',spacersl*725.0,
' km'
1342 print *,
'Unit of wavelength: Angstrom (0.1 nm) '
1344 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
1346 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
1348 write(*,
'(a,f8.1,a)')
' Supposed width of slit: ',wslit*725.0,
' km'
1353 print *,
'Unit of wavelength: Angstrom (0.1 nm) '
1355 write(*,
'(a,f5.3,a,f5.1,a)')
' Pixel: ',wlrsl,
' Angstrom x ',spacersl*725.0,
' km'
1356 print *,
'Unit of length: arcsec (~725 km)'
1357 write(*,
'(a,f8.1,a)')
' Location of slit: xI1 = ',
location_slit,
' arcsec'
1358 write(*,
'(a,f8.1,a)')
' Width of slit: ',wslit,
' arcsec'
1361 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
1363 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
1365 write(*,
'(a,f8.1,a)')
' Location of slit: xI1 = ',
location_slit,
' Unit_length'
1366 write(*,
'(a,f8.1,a)')
' Width of slit: ',wslit*725.d0,
' km'
1369 if (
mype==0) print *,
'Direction of the slit: parallel to xI2 vector'
1370 if (coordinate==cartesian .or. coordinate==spherical)
then
1373 call mpistop(
"EUV spectrum synthesis: support for sperical coordinates is to be added!")
1377 if (
mype==0) print *,
'###################################################'
1383 integer,
intent(in) :: qunit
1384 character(20),
intent(in) :: datatype
1387 integer :: numWL,numXS,iwL,ixS,numWI,numS
1388 double precision :: dwLg,xSmin,xSmax,wLmin,wLmax
1389 double precision,
allocatable :: wL(:),xS(:),dwL(:),dxS(:)
1390 double precision,
allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
1391 integer :: strtype,nstrb,nbb,nuni,nstr,bnx
1392 double precision :: qs,dxfirst,dxmid,lenstr
1394 integer :: iigrid,igrid,j,dir_loc
1395 double precision :: xbmin(1:ndim),xbmax(1:ndim)
1398 numwl=4*int((spectrum_window_max-spectrum_window_min)/(4.d0*dwlg))
1399 wlmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwlg*numwl/2
1400 wlmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwlg*numwl/2
1401 allocate(wl(numwl),dwl(numwl))
1404 wl(iwl)=wlmin+iwl*dwlg-half*dwlg
1407 select case(direction_slit)
1409 numxs=domain_nx1*2**(refine_max_level-1)
1414 strtype=stretch_type(1)
1415 nstrb=nstretchedblocks_baselevel(1)
1416 qs=qstretch_baselevel(1)
1417 if (mype==0) print *,
'Direction of the slit: x'
1419 numxs=domain_nx2*2**(refine_max_level-1)
1424 strtype=stretch_type(2)
1425 nstrb=nstretchedblocks_baselevel(2)
1426 qs=qstretch_baselevel(2)
1427 if (mype==0) print *,
'Direction of the slit: y'
1429 numxs=domain_nx3*2**(refine_max_level-1)
1434 strtype=stretch_type(3)
1435 nstrb=nstretchedblocks_baselevel(3)
1436 qs=qstretch_baselevel(3)
1437 if (mype==0) print *,
'Direction of the slit: z'
1439 call mpistop(
'Wrong direction_slit')
1442 allocate(xs(numxs),dxs(numxs),spectra(numwl,numxs),spectra_rc(numwl,numxs))
1444 allocate(wi(numwl,numxs,numwi))
1446 select case(strtype)
1448 dxs(:)=(xsmax-xsmin)/numxs
1450 xs(ixs)=xsmin+dxs(ixs)*(ixs-half)
1453 qs=qs**(one/2**(refine_max_level-1))
1454 dxfirst=(xsmax-xsmin)*(one-qs)/(one-qs**numxs)
1457 dxs(ixs)=dxfirst*qs**(ixs-1)
1458 xs(ixs)=dxs(1)/(one-qs)*(one-qs**(ixs-1))+half*dxs(ixs)
1464 lenstr=(xsmax-xsmin)/(2.d0+nuni*(one-qs)/(one-qs**nstr))
1465 dxfirst=(xsmax-xsmin)/(dble(nuni)+2.d0/(one-qs)*(one-qs**nstr))
1468 nstr=nstr*2**(refine_max_level-1)
1469 nuni=nuni*2**(refine_max_level-1)
1470 qs=qs**(one/2**(refine_max_level-1))
1471 dxfirst=lenstr*(one-qs)/(one-qs**nstr)
1472 dxmid=dxmid/2**(refine_max_level-1)
1474 if(nuni .gt. 0)
then
1475 do ixs=nstr+1,nstr+nuni
1477 xs(ixs)=lenstr+(dble(ixs)-0.5d0-nstr)*dxs(ixs)+xsmin
1482 dxs(ixs)=dxfirst*qs**(nstr-ixs)
1483 xs(ixs)=xsmin+lenstr-dxs(ixs)*half-dxfirst*(one-qs**(nstr-ixs))/(one-qs)
1486 do ixs=nstr+nuni+1,numxs
1487 dxs(ixs)=dxfirst*qs**(ixs-nstr-nuni-1)
1488 xs(ixs)=xsmax-lenstr+dxs(ixs)*half+dxfirst*(one-qs**(ixs-nstr-nuni-1))/(one-qs)
1491 call mpistop(
"unknown stretch type")
1494 if (los_phi==0 .and. los_theta==90 .and. direction_slit==2)
then
1497 else if (los_phi==0 .and. los_theta==90 .and. direction_slit==3)
then
1500 else if (los_phi==90 .and. los_theta==90 .and. direction_slit==1)
then
1503 else if (los_phi==90 .and. los_theta==90 .and. direction_slit==3)
then
1506 else if (los_theta==0 .and. direction_slit==1)
then
1509 else if (los_theta==0 .and. direction_slit==2)
then
1513 call mpistop(
'Wrong combination of LOS and slit direction!')
1516 if (dir_loc==1)
then
1517 if (location_slit>xprobmax1 .or. location_slit<xprobmin1)
then
1518 call mpistop(
'Wrong value for location_slit!')
1520 if(mype==0)
write(*,
'(a,f8.1,a)')
' Location of slit: x = ',location_slit,
' Unit_length'
1521 else if (dir_loc==2)
then
1522 if (location_slit>xprobmax2 .or. location_slit<xprobmin2)
then
1523 call mpistop(
'Wrong value for location_slit!')
1525 if(mype==0)
write(*,
'(a,f8.1,a)')
' Location of slit: y = ',location_slit,
' Unit_length'
1527 if (location_slit>xprobmax3 .or. location_slit<xprobmin3)
then
1528 call mpistop(
'Wrong value for location_slit!')
1530 if(mype==0)
write(*,
'(a,f8.1,a)')
' Location of slit: z = ',location_slit,
' Unit_length'
1535 do iigrid=1,igridstail; igrid=igrids(iigrid);
1536 ^d&xbmin(^d)=rnode(rpxmin^d_,igrid);
1537 ^d&xbmax(^d)=rnode(rpxmax^d_,igrid);
1538 if (location_slit>=xbmin(dir_loc) .and. location_slit<xbmax(dir_loc))
then
1544 call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
1545 mpi_sum,icomm,ierrmpi)
1548 if (spectra_rc(iwl,ixs)>smalldouble)
then
1549 wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
1556 call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
1558 deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
1565 integer,
intent(in) :: igrid,numWL,numXS,dir_loc
1567 double precision,
intent(in) :: wL(numWL),dwL(numWL)
1568 double precision,
intent(inout) :: spectra(numWL,numXS)
1570 integer :: direction_LOS
1571 integer :: ixO^L,ixI^L,ix^D,ixOnew
1572 double precision,
allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
1573 double precision :: wlc,wlwd
1576 double precision :: logTe,lineCent
1577 character (30) :: ion
1578 double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
1580 integer :: levelg,rft,ixSmin,ixSmax,iwL
1581 double precision :: flux_pix,dL
1583 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1585 if (los_phi==0 .and. los_theta==90)
then
1587 else if (los_phi==90 .and. los_theta==90)
then
1593 ^d&ixomin^d=ixmlo^d\
1594 ^d&ixomax^d=ixmhi^d\
1595 ^d&iximin^d=ixglo^d\
1596 ^d&iximax^d=ixghi^d\
1597 allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
1600 if (dir_loc==1)
then
1601 do ix1=ixomin1,ixomax1
1602 if (location_slit>=(ps(igrid)%x(ix^d,1)-
half*ps(igrid)%dx(ix^d,1)) .and. &
1603 location_slit<(ps(igrid)%x(ix^d,1)+
half*ps(igrid)%dx(ix^d,1)))
then
1609 else if (dir_loc==2)
then
1610 do ix2=ixomin2,ixomax2
1611 if (location_slit>=(ps(igrid)%x(ix^d,2)-
half*ps(igrid)%dx(ix^d,2)) .and. &
1612 location_slit<(ps(igrid)%x(ix^d,2)+
half*ps(igrid)%dx(ix^d,2)))
then
1619 do ix3=ixomin3,ixomax3
1620 if (location_slit>=(ps(igrid)%x(ix^d,3)-
half*ps(igrid)%dx(ix^d,3)) .and. &
1621 location_slit<(ps(igrid)%x(ix^d,3)+
half*ps(igrid)%dx(ix^d,3)))
then
1629 call get_euv(spectrum_wl,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
1630 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
1631 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
1632 v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
1633 call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
1634 call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
1635 te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
1638 levelg=ps(igrid)%level
1639 rft=2**(refine_max_level-levelg)
1641 {
do ix^d=ixomin^d,ixomax^d\}
1644 wlc=linecent*(1.d0+v(ix^d)*unit_velocity*1.d2/
const_c)
1646 wlc=linecent*(1.d0+v(ix^d)*unit_velocity/
const_c)
1648 wlwd=sqrt(
kb_cgs*te(ix^d)*unit_temperature/(mass*
mp_cgs))
1651 select case(direction_slit)
1653 ixsmin=(block_nx1*(node(pig1_,igrid)-1)+(ix1-ixomin1))*rft+1
1654 ixsmax=(block_nx1*(node(pig1_,igrid)-1)+(ix1-ixomin1+1))*rft
1656 ixsmin=(block_nx2*(node(pig2_,igrid)-1)+(ix2-ixomin2))*rft+1
1657 ixsmax=(block_nx2*(node(pig2_,igrid)-1)+(ix2-ixomin2+1))*rft
1659 ixsmin=(block_nx3*(node(pig3_,igrid)-1)+(ix3-ixomin3))*rft+1
1660 ixsmax=(block_nx3*(node(pig3_,igrid)-1)+(ix3-ixomin3+1))*rft
1663 select case(direction_los)
1665 dl=ps(igrid)%dx(ix^d,1)*unit_length
1667 dl=ps(igrid)%dx(ix^d,2)*unit_length
1669 dl=ps(igrid)%dx(ix^d,3)*unit_length
1671 if (si_unit) dl=dl*1.d2
1674 flux_pix=flux(ix^d)*wlrsl*dl*
exp_clamped(-(wl(iwl)-wlc)**2/(2*wlwd**2))/(sqrt(2*
dpi)*wlwd)
1676 flux_pix=flux_pix*wslit/spacersl
1677 spectra(iwl,ixsmin:ixsmax)=spectra(iwl,ixsmin:ixsmax)+flux_pix
1683 deallocate(flux,v,pth,te,rho)
1689 integer,
intent(in) :: qunit
1690 character(20),
intent(in) :: datatype
1693 integer :: numWL,numXS,iwL,ixS,numWI,ix^D
1694 double precision :: dwLg,dxSg,xSmin,xSmax,xScent,wLmin,wLmax
1695 double precision,
allocatable :: wL(:),xS(:),dwL(:),dxS(:)
1696 double precision,
allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
1697 double precision :: vec_cor(1:3),xI_cor(1:2)
1698 double precision :: res,r_loc,r_max
1701 character (30) :: ion
1702 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1703 double precision :: unitv,arcsec,RHESSI_rsl,pixel
1704 integer :: iigrid,igrid,i,j,numS
1705 double precision :: xLmin,xLmax,xslit
1707 if (coordinate==spherical)
then
1715 if (coordinate==spherical)
then
1716 xsmin=-abs(xprobmax1)
1717 xsmax=abs(xprobmax1)
1720 if (ix1==1) vec_cor(1)=xprobmin1
1721 if (ix1==2) vec_cor(1)=xprobmax1
1723 if (ix2==1) vec_cor(2)=xprobmin2
1724 if (ix2==2) vec_cor(2)=xprobmax2
1726 if (ix3==1) vec_cor(3)=xprobmin3
1727 if (ix3==2) vec_cor(3)=xprobmax3
1729 r_loc=(vec_cor(1)-x_origin(1))**2
1730 r_loc=r_loc+(vec_cor(2)-x_origin(2))**2
1731 r_loc=r_loc+(vec_cor(3)-x_origin(3))**2
1733 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1736 r_max=max(r_max,r_loc)
1740 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1744 xsmin=min(xsmin,xi_cor(2))
1745 xsmax=max(xsmax,xi_cor(2))
1756 xscent=(xsmin+xsmax)/2.d0
1760 arcsec=7.25d5/unit_length
1762 arcsec=7.25d7/unit_length
1764 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1765 dxsg=spacersl*arcsec
1766 numxs=ceiling((xsmax-xscent)/dxsg)
1767 xsmin=xscent-numxs*dxsg
1768 xsmax=xscent+numxs*dxsg
1771 numwl=2*int((spectrum_window_max-spectrum_window_min)/(2.d0*dwlg))
1772 wlmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwlg*numwl/2
1773 wlmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwlg*numwl/2
1774 allocate(wl(numwl),dwl(numwl),xs(numxs),dxs(numxs))
1776 allocate(wi(numwl,numxs,numwi),spectra(numwl,numxs),spectra_rc(numwl,numxs))
1778 wl(iwl)=wlmin+iwl*dwlg-half*dwlg
1782 xs(ixs)=xsmin+dxsg*(ixs-half)
1788 do iigrid=1,igridstail; igrid=igrids(iigrid);
1790 if (ix1==1) vec_cor(1)=rnode(rpxmin1_,igrid)
1791 if (ix1==2) vec_cor(1)=rnode(rpxmax1_,igrid)
1793 if (ix2==1) vec_cor(2)=rnode(rpxmin2_,igrid)
1794 if (ix2==2) vec_cor(2)=rnode(rpxmax2_,igrid)
1796 if (ix3==1) vec_cor(3)=rnode(rpxmin3_,igrid)
1797 if (ix3==2) vec_cor(3)=rnode(rpxmax3_,igrid)
1799 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1803 xlmin=min(xlmin,xi_cor(1))
1804 xlmax=max(xlmax,xi_cor(1))
1810 if (activate_unit_arcsec)
then
1811 xslit=location_slit*arcsec
1815 if (xslit>=xlmin-wslit*arcsec .and. xslit<=xlmax+wslit*arcsec)
then
1821 call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
1822 mpi_sum,icomm,ierrmpi)
1825 if (spectra_rc(iwl,ixs)>smalldouble)
then
1826 wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
1833 if (activate_unit_arcsec)
then
1838 call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
1840 deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
1846 integer,
intent(in) :: igrid,numWL,numXS
1847 double precision,
intent(in) :: wL(numWL),xS(numXS)
1848 double precision,
intent(in) :: dwLg,dxSg
1849 double precision,
intent(inout) :: spectra(numWL,numXS)
1852 integer :: ixO^L,ixI^L,ix^D,ixOnew,j
1853 double precision,
allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
1854 double precision :: wlc,wlwd,res,dst_slit,xslit,arcsec
1855 double precision :: vloc(1:3),xloc(1:3),dxloc(1:3),xIloc(1:2),dxIloc(1:2)
1856 integer :: nSubC^D,iSubC^D,iwL,ixS,ixSmin,ixSmax,iwLmin,iwLmax,nwL
1857 double precision :: slit_width,dxSubC^D,xerf^L,fluxSubC
1858 double precision :: xSubC(1:3),xCent(1:2)
1861 double precision :: logTe,lineCent
1862 character (30) :: ion
1863 double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
1864 double precision :: sigma_wl,sigma_xs,factor
1867 arcsec=7.25d5/unit_length
1869 arcsec=7.25d7/unit_length
1871 if (activate_unit_arcsec)
then
1872 xslit=location_slit*arcsec
1877 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1879 ^d&ixomin^d=ixmlo^d\
1880 ^d&ixomax^d=ixmhi^d\
1881 ^d&iximin^d=ixglo^d\
1882 ^d&iximax^d=ixghi^d\
1883 allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
1885 call get_euv(spectrum_wl,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
1886 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
1887 call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
1888 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
1889 call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
1890 te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
1891 {
do ix^d=ixomin^d,ixomax^d\}
1893 vloc(j)=ps(igrid)%w(ix^d,iw_mom(j))/rho(ix^d)
1901 slit_width=wslit*arcsec
1902 sigma_wl=sigma_psf*dwlg
1903 sigma_xs=sigma_psf*dxsg
1904 {
do ix^d=ixomin^d,ixomax^d\}
1905 if (flux(ix^d)>smalldouble)
then
1906 xloc(1:3)=ps(igrid)%x(ix^d,1:3)
1907 dxloc(1:3)=ps(igrid)%dx(ix^d,1:3)
1911 if (xiloc(1)>=xslit-half*(slit_width+dxiloc(1)) .and. &
1912 xiloc(1)<=xslit+half*(slit_width+dxiloc(1)))
then
1914 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi1(^d))/(slit_width/16.d0)));
1915 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi2(^d))/(dxsg/4.d0)));
1916 ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
1919 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length*1.d2/dxsg/dxsg
1920 wlc=linecent*(1.d0+v(ix^d)*unit_velocity*1.d2/const_c)
1922 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length/dxsg/dxsg
1923 wlc=linecent*(1.d0+v(ix^d)*unit_velocity/const_c)
1925 wlwd=sqrt(kb_cgs*te(ix^d)*unit_temperature/(mass*mp_cgs))
1926 wlwd=wlwd*linecent/const_c
1928 {
do isubc^d=1,nsubc^d\}
1929 ^d&xsubc(^d)=xloc(^d)-half*dxloc(^d)+(isubc^d-half)*dxsubc^d;
1931 dst_slit=abs(xcent(1)-xslit)
1932 if (dst_slit<=half*slit_width)
then
1933 ixs=floor((xcent(2)-(xs(1)-half*dxsg))/dxsg)+1
1935 ixsmax=min(ixs+3,numxs)
1936 iwl=floor((wlc-(wl(1)-half*dwlg))/dwlg)+1
1937 nwl=3*ceiling(wlwd/dwlg+1)
1938 iwlmin=max(1,iwl-nwl)
1939 iwlmax=min(iwl+nwl,numwl)
1941 do iwl=iwlmin,iwlmax
1942 do ixs=ixsmin,ixsmax
1943 xerfmin1=(wl(iwl)-half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
1944 xerfmax1=(wl(iwl)+half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
1945 xerfmin2=(xs(ixs)-half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
1946 xerfmax2=(xs(ixs)+half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
1947 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
1948 spectra(iwl,ixs)=spectra(iwl,ixs)+fluxsubc*factor
1958 deallocate(flux,v,pth,te)
1966 integer,
intent(in) :: qunit
1968 character(20) :: datatype
1971 character (30) :: ion
1972 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1973 double precision :: t0,t1
1976 datatype=
'image_euv'
1981 print *,
'###################################################'
1982 print *,
'Systhesizing EUV image'
1983 write(*,
'(a,f8.3,a)')
' Wavelength: ',linecent,
' Angstrom'
1984 print *,
'Unit of EUV flux: DN s^-1 pixel^-1'
1989 call mpistop(
'EUV dat-resolution needs Cartesian or spherical native rays')
1991 write(*,
'(a,f7.1,a,f7.1,a,f5.1,a,f5.1,a)')
' Supposed Pixel: ',spacersl*725.0,
' km x ',spacersl*725.0, &
1992 ' km (', spacersl,
' arcsec x ', spacersl,
' arcsec)'
1994 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
1996 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2010 call mpistop(
'ERROR: Wrong LOS for synthesizing emission!')
2014 write(*,
'(a,f7.1,a,f7.1,a,f5.1,a,f5.1,a)')
' Pixel: ',spacersl*725.0,
' km x ',spacersl*725.0,
' km (', &
2015 spacersl,
' arcsec x ', spacersl,
' arcsec)'
2017 print *,
'Unit of length: arcsec (~725 km)'
2020 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2022 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2026 if (coordinate==cartesian)
then
2028 '] of the simulation box is located at [X=0,Y=0] of the image'
2030 else if (coordinate==spherical)
then
2031 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'
2034 call mpistop(
"EUV synthesis: this coordinate is not supported!")
2039 if (
mype==0) print *,
'time comsuming: ',t1-t0,
' s'
2040 if (
mype==0) print *,
'###################################################'
2047 integer,
intent(in) :: qunit
2049 character(20) :: datatype
2050 double precision :: RHESSI_rsl
2051 double precision :: t0,t1
2054 datatype=
'image_sxr'
2059 print *,
'###################################################'
2060 print *,
'Systhesizing SXR image (observed at 1 AU).'
2065 if (coordinate/=cartesian)
call mpistop(
'SXR synthesis: only cartesian is supported for .dat resolution!')
2067 print *,
'Unit of SXR flux: photons cm^-2 s^-1 pixel^-1'
2068 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, &
2069 ' Mm (', rhessi_rsl,
' arcsec x ', rhessi_rsl,
' arcsec)'
2071 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2073 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2083 call mpistop(
'ERROR: Wrong LOS for synthesizing emission!')
2087 print *,
'Unit of SXR flux: photons cm^-2 s^-1 pixel^-1'
2088 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, &
2089 ' Mm (', rhessi_rsl,
' arcsec x ', rhessi_rsl,
' arcsec)'
2091 print *,
'Unit of length: arcsec (~725 km)'
2094 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2096 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2100 if (coordinate==cartesian)
then
2102 '] of the simulation box is located at [X=0,Y=0] of the image'
2104 else if (coordinate==spherical)
then
2105 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'
2108 call mpistop(
"SXR synthesis: this coordinate is not supported!")
2113 if (
mype==0) print *,
'time comsuming:',t1-t0
2114 if (
mype==0) print *,
'###################################################'
2121 integer,
intent(in) :: qunit
2123 character(20) :: datatype
2124 double precision :: LASCO_rsl
2126 if (
mype==0) print *,
'###################################################'
2130 if (
mype==0) print *,
'Systhesizing white light image (observed by LASCO/C1).'
2133 if (
mype==0) print *,
'Systhesizing white light image (observed by LASCO/C2).'
2136 if (
mype==0) print *,
'Systhesizing white light image (observed by LASCO/C3).'
2138 call mpistop(
'Whitelight synthesis: instrument is not supported!')
2141 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 (', &
2142 lasco_rsl,
' arcsec x ', lasco_rsl,
' arcsec) '
2143 if (
mype==0) print *,
'Unit of white light flux: average Sun brightness'
2145 datatype=
'image_whitelight'
2150 print *,
'Unit of length: arcsec (~725 km)'
2153 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2155 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2160 if (coordinate==spherical)
then
2161 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'
2164 call mpistop(
"Whitelight synthesis: this coordinate is not supported!")
2167 if (
mype==0) print *,
'###################################################'
2172 EUV,Dpl,nOut1,nOut2,xOut1,xOut2,&
2173 dxOut1,dxOut2,wOut,numWOut,Tau,EUVthin)
2177 integer,
intent(in) :: nSrc1,nSrc2
2178 double precision,
intent(in) :: xSrc1(nSrc1),xSrc2(nSrc2)
2179 double precision,
intent(in) :: dxSrc1(nSrc1),dxSrc2(nSrc2)
2180 double precision,
intent(in) :: EUV(nSrc1,nSrc2),Dpl(nSrc1,nSrc2)
2181 integer,
intent(out) :: nOut1,nOut2,numWOut
2182 double precision,
allocatable,
intent(out) :: xOut1(:),xOut2(:),dxOut1(:),dxOut2(:)
2183 double precision,
allocatable,
intent(out) :: wOut(:,:,:)
2184 double precision,
intent(in),
optional :: Tau(nSrc1,nSrc2),EUVthin(nSrc1,nSrc2)
2186 integer :: mass,ixS1,ixS2,ixP1,ixP2,ixC1,ixC2,iw
2187 integer :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
2188 character(30) :: ion
2189 double precision :: logTe,lineCent,spaceRsl,wlRsl,sigma_PSF,wslit
2190 double precision :: arcsec,dxInst,xMin1,xMax1,xMin2,xMax2,xCent1,xCent2
2191 double precision :: sigma0,xerfmin1,xerfmax1,xerfmin2,xerfmax2
2192 double precision :: factor,weightSum,weightNorm,thinVal,tauVal
2193 double precision,
allocatable :: dplNum(:,:),thinOut(:,:),tauOut(:,:),tauWeight(:,:)
2201 dxinst=spacersl*arcsec
2202 if (dxinst<=
zero)
call mpistop(
"instrument_postprocess has non-positive pixel size")
2204 xmin1=minval(xsrc1-
half*dxsrc1)
2205 xmax1=maxval(xsrc1+
half*dxsrc1)
2206 xmin2=minval(xsrc2-
half*dxsrc2)
2207 xmax2=maxval(xsrc2+
half*dxsrc2)
2208 xcent1=
half*(xmin1+xmax1)
2209 xcent2=
half*(xmin2+xmax2)
2210 nout1=16*max(1,ceiling((xmax1-xmin1)/(16.d0*dxinst)))
2211 nout2=16*max(1,ceiling((xmax2-xmin2)/(16.d0*dxinst)))
2212 xmin1=xcent1-
half*dble(nout1)*dxinst
2213 xmin2=xcent2-
half*dble(nout2)*dxinst
2215 allocate(xout1(nout1),xout2(nout2),dxout1(nout1),dxout2(nout2))
2217 xout1(ixp1)=xmin1+dxinst*(dble(ixp1)-
half)
2221 xout2(ixp2)=xmin2+dxinst*(dble(ixp2)-
half)
2226 if (
present(tau) .and.
output_tau) numwout=numwout+1
2228 allocate(wout(nout1,nout2,numwout),dplnum(nout1,nout2))
2231 if (
present(euvthin))
then
2232 allocate(thinout(nout1,nout2))
2236 allocate(tauout(nout1,nout2),tauweight(nout1,nout2))
2241 sigma0=sigma_psf*dxinst
2246 if (
present(euvthin)) thinval=euvthin(ixs1,ixs2)
2247 if (
present(tau)) tauval=tau(ixs1,ixs2)
2251 ixc1=floor((xsrc1(ixs1)-(xout1(1)-
half*dxinst))/dxinst)+1
2252 ixc2=floor((xsrc2(ixs2)-(xout2(1)-
half*dxinst))/dxinst)+1
2253 ixpmin1=max(1,ixc1-3)
2254 ixpmax1=min(nout1,ixc1+3)
2255 ixpmin2=max(1,ixc2-3)
2256 ixpmax2=min(nout2,ixc2+3)
2259 do ixp1=ixpmin1,ixpmax1
2260 do ixp2=ixpmin2,ixpmax2
2261 xerfmin1=((xout1(ixp1)-
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2262 xerfmax1=((xout1(ixp1)+
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2263 xerfmin2=((xout2(ixp2)-
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2264 xerfmax2=((xout2(ixp2)+
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2265 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
2266 weightsum=weightsum+factor
2269 if (weightsum<=
zero) cycle
2271 do ixp1=ixpmin1,ixpmax1
2272 do ixp2=ixpmin2,ixpmax2
2273 xerfmin1=((xout1(ixp1)-
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2274 xerfmax1=((xout1(ixp1)+
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2275 xerfmin2=((xout2(ixp2)-
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2276 xerfmax2=((xout2(ixp2)+
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2277 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
2278 weightnorm=factor/weightsum
2279 wout(ixp1,ixp2,1)=wout(ixp1,ixp2,1)+euv(ixs1,ixs2)*weightnorm
2280 dplnum(ixp1,ixp2)=dplnum(ixp1,ixp2)+euv(ixs1,ixs2)*dpl(ixs1,ixs2)*weightnorm
2281 if (
present(euvthin)) thinout(ixp1,ixp2)=thinout(ixp1,ixp2)+thinval*weightnorm
2283 tauout(ixp1,ixp2)=tauout(ixp1,ixp2)+tauval*weightnorm
2284 tauweight(ixp1,ixp2)=tauweight(ixp1,ixp2)+weightnorm
2294 wout(ixp1,ixp2,2)=dplnum(ixp1,ixp2)/wout(ixp1,ixp2,1)
2296 wout(ixp1,ixp2,2)=
zero
2305 if (tauweight(ixp1,ixp2)>
zero)
then
2306 wout(ixp1,ixp2,iw)=tauout(ixp1,ixp2)/tauweight(ixp1,ixp2)
2308 wout(ixp1,ixp2,iw)=
zero
2318 wout(ixp1,ixp2,iw)=min(
one,max(
zero,(thinout(ixp1,ixp2)-wout(ixp1,ixp2,1))/thinout(ixp1,ixp2)))
2320 wout(ixp1,ixp2,iw)=
zero
2327 write(*,
'(a,2(i8,1x),a,2(i8,1x),a,1pe12.5)') &
2328 ' instrument_postprocess EUV grid src/out: ',nsrc1,nsrc2,
' -> ',nout1,nout2,
' dx=',dxinst
2332 if (
allocated(thinout))
deallocate(thinout)
2333 if (
allocated(tauout))
deallocate(tauout,tauweight)
2337 Bright,nOut1,nOut2,xOut1,xOut2,&
2338 dxOut1,dxOut2,wOut,numWOut,Tau,BrightThin)
2341 integer,
intent(in) :: nSrc1,nSrc2
2342 double precision,
intent(in) :: xSrc1(nSrc1),xSrc2(nSrc2)
2343 double precision,
intent(in) :: dxSrc1(nSrc1),dxSrc2(nSrc2)
2344 double precision,
intent(in) :: Bright(nSrc1,nSrc2)
2345 integer,
intent(out) :: nOut1,nOut2,numWOut
2346 double precision,
allocatable,
intent(out) :: xOut1(:),xOut2(:),dxOut1(:),dxOut2(:)
2347 double precision,
allocatable,
intent(out) :: wOut(:,:,:)
2348 double precision,
intent(in),
optional :: Tau(nSrc1,nSrc2),BrightThin(nSrc1,nSrc2)
2350 integer :: ixS1,ixS2,ixP1,ixP2,ixC1,ixC2,iw,nStencil
2351 integer :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
2352 double precision :: arcsec,beamPixel,beamSigma,xMin1,xMax1,xMin2,xMax2,xCent1,xCent2
2353 double precision :: distance1,distance2,weight,cellArea,thinVal,tauVal
2354 double precision,
allocatable :: norm(:,:),thinOut(:,:),tauOut(:,:),tauNorm(:,:)
2367 if (beamsigma<=zero .or. beampixel<=zero)
then
2368 call mpistop(
"radio beam postprocess has non-positive beam or pixel size")
2371 xmin1=minval(xsrc1-half*dxsrc1)
2372 xmax1=maxval(xsrc1+half*dxsrc1)
2373 xmin2=minval(xsrc2-half*dxsrc2)
2374 xmax2=maxval(xsrc2+half*dxsrc2)
2375 xcent1=half*(xmin1+xmax1)
2376 xcent2=half*(xmin2+xmax2)
2377 nout1=16*max(1,ceiling((xmax1-xmin1)/(16.d0*beampixel)))
2378 nout2=16*max(1,ceiling((xmax2-xmin2)/(16.d0*beampixel)))
2379 xmin1=xcent1-half*dble(nout1)*beampixel
2380 xmin2=xcent2-half*dble(nout2)*beampixel
2382 allocate(xout1(nout1),xout2(nout2),dxout1(nout1),dxout2(nout2))
2384 xout1(ixp1)=xmin1+beampixel*(dble(ixp1)-half)
2385 dxout1(ixp1)=beampixel
2388 xout2(ixp2)=xmin2+beampixel*(dble(ixp2)-half)
2389 dxout2(ixp2)=beampixel
2393 if (
present(tau) .and.
output_tau) numwout=numwout+1
2395 allocate(wout(nout1,nout2,numwout),norm(nout1,nout2))
2398 if (
present(brightthin))
then
2399 allocate(thinout(nout1,nout2))
2403 allocate(tauout(nout1,nout2),taunorm(nout1,nout2))
2408 nstencil=max(3,ceiling(4.d0*beamsigma/beampixel)+1)
2413 if (
present(brightthin)) thinval=brightthin(ixs1,ixs2)
2414 if (
present(tau)) tauval=tau(ixs1,ixs2)
2415 if (abs(bright(ixs1,ixs2))<=smalldouble .and. abs(thinval)<=smalldouble .and. &
2416 abs(tauval)<=smalldouble) cycle
2418 ixc1=floor((xsrc1(ixs1)-(xout1(1)-half*beampixel))/beampixel)+1
2419 ixc2=floor((xsrc2(ixs2)-(xout2(1)-half*beampixel))/beampixel)+1
2420 ixpmin1=max(1,ixc1-nstencil)
2421 ixpmax1=min(nout1,ixc1+nstencil)
2422 ixpmin2=max(1,ixc2-nstencil)
2423 ixpmax2=min(nout2,ixc2+nstencil)
2424 cellarea=max(smalldouble,dxsrc1(ixs1)*dxsrc2(ixs2))
2426 do ixp1=ixpmin1,ixpmax1
2427 distance1=xout1(ixp1)-xsrc1(ixs1)
2428 do ixp2=ixpmin2,ixpmax2
2429 distance2=xout2(ixp2)-xsrc2(ixs2)
2430 weight=
exp_clamped(-half*(distance1**2+distance2**2)/beamsigma**2)*cellarea
2431 wout(ixp1,ixp2,1)=wout(ixp1,ixp2,1)+bright(ixs1,ixs2)*weight
2432 norm(ixp1,ixp2)=norm(ixp1,ixp2)+weight
2433 if (
present(brightthin)) thinout(ixp1,ixp2)=thinout(ixp1,ixp2)+thinval*weight
2435 tauout(ixp1,ixp2)=tauout(ixp1,ixp2)+tauval*weight
2436 taunorm(ixp1,ixp2)=taunorm(ixp1,ixp2)+weight
2445 if (norm(ixp1,ixp2)>zero)
then
2446 wout(ixp1,ixp2,1)=wout(ixp1,ixp2,1)/norm(ixp1,ixp2)
2447 if (
present(brightthin)) thinout(ixp1,ixp2)=thinout(ixp1,ixp2)/norm(ixp1,ixp2)
2449 wout(ixp1,ixp2,1)=zero
2450 if (
present(brightthin)) thinout(ixp1,ixp2)=zero
2460 if (taunorm(ixp1,ixp2)>zero)
then
2461 wout(ixp1,ixp2,iw)=tauout(ixp1,ixp2)/taunorm(ixp1,ixp2)
2463 wout(ixp1,ixp2,iw)=zero
2472 if (thinout(ixp1,ixp2)>smalldouble)
then
2473 wout(ixp1,ixp2,iw)=min(one,max(zero,(thinout(ixp1,ixp2)-wout(ixp1,ixp2,1))/thinout(ixp1,ixp2)))
2475 wout(ixp1,ixp2,iw)=zero
2482 write(*,
'(a,2(i8,1x),a,2(i8,1x),a,2(1pe12.5,1x))') &
2483 ' radio_beam_postprocess grid src/out: ',nsrc1,nsrc2,
' -> ',nout1,nout2,&
2488 if (
allocated(thinout))
deallocate(thinout)
2489 if (
allocated(tauout))
deallocate(tauout,taunorm)
2498 integer,
intent(in) :: qunit
2499 character(20),
intent(in) :: datatype
2502 double precision :: dx^D
2503 integer :: numX^D,ix^D
2504 double precision,
allocatable :: EUV(:,:),EUVs(:,:),Dpl(:,:),Dpls(:,:)
2505 double precision,
allocatable :: EUVthin(:,:),Tau(:,:)
2506 double precision,
allocatable :: SXR(:,:),SXRs(:,:),wI(:,:,:)
2507 double precision,
allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:),dxIi
2508 integer :: numXI1,numXI2,numSI,numWI,iw
2509 double precision :: xI^L
2510 integer :: iigrid,igrid,i,j
2511 double precision,
allocatable :: xIF1(:),xIF2(:),dxIF1(:),dxIF2(:)
2512 double precision,
allocatable :: xIP1(:),xIP2(:),dxIP1(:),dxIP2(:),wIP(:,:,:)
2513 integer :: nXIF1,nXIF2
2514 integer :: nXIP1,nXIP2,numWIP
2515 double precision :: xIF^L
2516 double precision :: vec_cor(1:3),xI_cor(1:2),dxDDA,xIcent1,xIcent2
2518 double precision :: unitv,arcsec,RHESSI_rsl
2519 integer :: strtype^D,nstrb^D,nbb^D,nuni^D,nstr^D,bnx^D
2520 double precision :: qs^D,dxfirst^D,dxmid^D,lenstr^D
2521 logical :: has_doppler_output,has_thick_output
2534 xicent1=
half*(xifmin1+xifmax1)
2535 xicent2=
half*(xifmin2+xifmax2)
2536 nxif1=max(1,ceiling((xifmax1-xifmin1)/dxdda))
2537 nxif2=max(1,ceiling((xifmax2-xifmin2)/dxdda))
2538 xifmin1=xicent1-
half*dble(nxif1)*dxdda
2539 xifmax1=xicent1+
half*dble(nxif1)*dxdda
2540 xifmin2=xicent2-
half*dble(nxif2)*dxdda
2541 xifmax2=xicent2+
half*dble(nxif2)*dxdda
2552 if (
mype==0)
write(*,
'(a,1pe12.5,a,2(i8,1x))') &
2553 ' spherical native dat-resolution image-plane dx: ',dxdda,
' n=',nxif1,nxif2
2558 if (ix1==1) vec_cor(1)=xprobmin1
2559 if (ix1==2) vec_cor(1)=xprobmax1
2561 if (ix2==1) vec_cor(2)=xprobmin2
2562 if (ix2==2) vec_cor(2)=xprobmax2
2564 if (ix3==1) vec_cor(3)=xprobmin3
2565 if (ix3==2) vec_cor(3)=xprobmax3
2567 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
2573 xifmin1=min(xifmin1,xi_cor(1))
2574 xifmax1=max(xifmax1,xi_cor(1))
2575 xifmin2=min(xifmin2,xi_cor(2))
2576 xifmax2=max(xifmax2,xi_cor(2))
2581 dxdda=min((xprobmax1-xprobmin1)/dble(numx1),&
2582 (xprobmax2-xprobmin2)/dble(numx2),&
2583 (xprobmax3-xprobmin3)/dble(numx3))
2584 xicent1=
half*(xifmin1+xifmax1)
2585 xicent2=
half*(xifmin2+xifmax2)
2586 nxif1=max(1,ceiling((xifmax1-xifmin1)/dxdda))
2587 nxif2=max(1,ceiling((xifmax2-xifmin2)/dxdda))
2588 xifmin1=xicent1-
half*dble(nxif1)*dxdda
2589 xifmax1=xicent1+
half*dble(nxif1)*dxdda
2590 xifmin2=xicent2-
half*dble(nxif2)*dxdda
2591 xifmax2=xicent2+
half*dble(nxif2)*dxdda
2619 if (
mype==0)
write(*,
'(a)')
' LOS vector: [-1.00 0.00 0.00]'
2620 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [ 0.00 1.00 0.00]'
2621 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 0.00 1.00]'
2639 if (
mype==0)
write(*,
'(a)')
' LOS vector: [ 0.00 -1.00 0.00]'
2640 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [-1.00 0.00 0.00]'
2641 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 0.00 1.00]'
2659 if (
mype==0)
write(*,
'(a)')
' LOS vector: [ 0.00 0.00 -1.00]'
2660 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [ 1.00 0.00 0.00]'
2661 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 1.00 0.00]'
2663 allocate(xif1(nxif1),xif2(nxif2),dxif1(nxif1),dxif2(nxif2))
2666 select case(strtype1)
2668 dxif1(:)=(xifmax1-xifmin1)/nxif1
2670 xif1(ix1)=xifmin1+dxif1(ix1)*(ix1-
half)
2674 dxfirst1=(xifmax1-xifmin1)*(
one-qs1)/(
one-qs1**nxif1)
2677 dxif1(ix1)=dxfirst1*qs1**(ix1-1)
2678 xif1(ix1)=dxif1(1)/(
one-qs1)*(
one-qs1**(ix1-1))+
half*dxif1(ix1)
2683 nuni1=nbb1-nstrb1*bnx1
2684 lenstr1=(xifmax1-xifmin1)/(2.d0+nuni1*(
one-qs1)/(
one-qs1**nstr1))
2685 dxfirst1=(xifmax1-xifmin1)/(dble(nuni1)+2.d0/(
one-qs1)*(
one-qs1**nstr1))
2691 dxfirst1=lenstr1*(
one-qs1)/(
one-qs1**nstr1)
2694 if(nuni1 .gt. 0)
then
2695 do ix1=nstr1+1,nstr1+nuni1
2697 xif1(ix1)=lenstr1+(dble(ix1)-0.5d0-nstr1)*dxif1(ix1)+xifmin1
2702 dxif1(ix1)=dxfirst1*qs1**(nstr1-ix1)
2703 xif1(ix1)=xifmin1+lenstr1-dxif1(ix1)*
half-dxfirst1*(
one-qs1**(nstr1-ix1))/(
one-qs1)
2706 do ix1=nstr1+nuni1+1,nxif1
2707 dxif1(ix1)=dxfirst1*qs1**(ix1-nstr1-nuni1-1)
2708 xif1(ix1)=xifmax1-lenstr1+dxif1(ix1)*
half+dxfirst1*(
one-qs1**(ix1-nstr1-nuni1-1))/(
one-qs1)
2711 call mpistop(
"unknown stretch type")
2714 select case(strtype2)
2716 dxif2(:)=(xifmax2-xifmin2)/nxif2
2718 xif2(ix2)=xifmin2+dxif2(ix2)*(ix2-
half)
2722 dxfirst2=(xifmax2-xifmin2)*(
one-qs2)/(
one-qs2**nxif2)
2725 dxif2(ix2)=dxfirst2*qs2**(ix2-1)
2726 xif2(ix2)=dxif2(1)/(
one-qs2)*(
one-qs2**(ix2-1))+
half*dxif2(ix2)
2731 nuni2=nbb2-nstrb2*bnx2
2732 lenstr2=(xifmax2-xifmin2)/(2.d0+nuni2*(
one-qs2)/(
one-qs2**nstr2))
2733 dxfirst2=(xifmax2-xifmin2)/(dble(nuni2)+2.d0/(
one-qs2)*(
one-qs2**nstr2))
2739 dxfirst2=lenstr2*(
one-qs2)/(
one-qs2**nstr2)
2742 if(nuni2 .gt. 0)
then
2743 do ix2=nstr2+1,nstr2+nuni2
2745 xif2(ix2)=lenstr2+(dble(ix2)-0.5d0-nstr2)*dxif2(ix2)+xifmin2
2750 dxif2(ix2)=dxfirst2*qs2**(nstr2-ix2)
2751 xif2(ix2)=xifmin2+lenstr2-dxif2(ix2)*
half-dxfirst2*(
one-qs2**(nstr2-ix2))/(
one-qs2)
2754 do ix2=nstr2+nuni2+1,nxif2
2755 dxif2(ix2)=dxfirst2*qs2**(ix2-nstr2-nuni2-1)
2756 xif2(ix2)=xifmax2-lenstr2+dxif2(ix2)*
half+dxfirst2*(
one-qs2**(ix2-nstr2-nuni2-1))/(
one-qs2)
2759 call mpistop(
"unknown stretch type")
2763 if (datatype==
'image_euv')
then
2772 allocate(wi(nxif1,nxif2,numwi))
2773 allocate(euv(nxif1,nxif2),dpl(nxif1,nxif2))
2775 allocate(euvthin(nxif1,nxif2),tau(nxif1,nxif2))
2784 if (has_doppler_output)
then
2789 allocate(euvs(nxif1,nxif2),dpls(nxif1,nxif2))
2798 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision, &
2803 do iigrid=1,igridstail; igrid=igrids(iigrid);
2807 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision, &
2809 call mpi_allreduce(dpls,dpl,numsi,mpi_double_precision, &
2812 if (has_doppler_output)
then
2816 deallocate(euvs,dpls)
2818 if (has_thick_output)
then
2819 if (has_doppler_output)
then
2821 has_thick_output,dpl=dpl,tau=tau,euvthin=euvthin)
2824 has_thick_output,tau=tau,euvthin=euvthin)
2826 else if (has_doppler_output)
then
2828 has_thick_output,dpl=dpl)
2838 euv,nxip1,nxip2,xip1,xip2,&
2839 dxip1,dxip2,wip,numwip,tau=tau,brightthin=euvthin)
2842 euv,nxip1,nxip2,xip1,xip2,&
2843 dxip1,dxip2,wip,numwip)
2847 euv,dpl,nxip1,nxip2,xip1,xip2,&
2848 dxip1,dxip2,wip,numwip,tau=tau,euvthin=euvthin)
2851 euv,dpl,nxip1,nxip2,xip1,xip2,&
2852 dxip1,dxip2,wip,numwip)
2854 call output_data(qunit,xip1,xip2,dxip1,dxip2,wip,nxip1,nxip2,numwip,datatype)
2855 deallocate(xip1,xip2,dxip1,dxip2,wip)
2857 call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
2860 deallocate(wi,euv,dpl,euvthin,tau)
2862 deallocate(wi,euv,dpl)
2867 if (datatype==
'image_sxr')
then
2875 allocate(wi(nxif1,nxif2,numwi))
2876 allocate(sxrs(nxif1,nxif2),sxr(nxif1,nxif2))
2879 do iigrid=1,igridstail; igrid=igrids(iigrid);
2883 call mpi_allreduce(sxrs,sxr,numsi,mpi_double_precision, &
2886 sxr=sxr*(rhessi_rsl*arcsec)**2
2894 call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
2895 deallocate(wi,sxr,sxrs)
2898 deallocate(xif1,xif2,dxif1,dxif2)
2905 integer,
intent(in) :: igrid,nXIF1,nXIF2
2906 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
2907 double precision,
intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
2909 double precision,
intent(out) :: SXR(nXIF1,nXIF2)
2911 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
2912 double precision :: xb^L,xd^D
2913 double precision,
allocatable :: flux(:^D&),opacity(:^D&)
2914 double precision,
allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
2915 double precision,
allocatable :: SXRg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
2916 integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
2917 double precision :: SXRt,xc^L,xg^L,r2,area_1AU
2918 integer :: ixP^L,ixP^D
2919 integer :: direction_LOS
2929 ^d&ixomin^d=ixmlo^d\
2930 ^d&ixomax^d=ixmhi^d\
2931 ^d&iximin^d=
ixglo^d\
2932 ^d&iximax^d=
ixghi^d\
2936 allocate(flux(ixi^s))
2937 allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
2938 dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
2939 dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
2940 dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
2945 levelg=ps(igrid)%level
2949 select case(direction_los)
2960 allocate(sxrg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
2966 select case(direction_los)
2968 do ix2=ixomin2,ixomax2
2969 ixgmin1=(ix2-1)*rft+1
2971 do ix3=ixomin3,ixomax3
2972 ixgmin2=(ix3-1)*rft+1
2975 do ix1=ixomin1,ixomax1
2978 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
2982 do ix3=ixomin3,ixomax3
2983 ixgmin1=(ix3-1)*rft+1
2985 do ix1=ixomin1,ixomax1
2986 ixgmin2=(ix1-1)*rft+1
2989 do ix2=ixomin2,ixomax2
2992 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
2996 do ix1=ixomin1,ixomax1
2997 ixgmin1=(ix1-1)*rft+1
2999 do ix2=ixomin2,ixomax2
3000 ixgmin2=(ix2-1)*rft+1
3003 do ix3=ixomin3,ixomax3
3006 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
3016 select case(direction_los)
3018 ixgmin1=(ixomin2-1)*rft+1
3020 ixgmin2=(ixomin3-1)*rft+1
3023 ixgmin1=(ixomin3-1)*rft+1
3025 ixgmin2=(ixomin1-1)*rft+1
3028 ixgmin1=(ixomin1-1)*rft+1
3030 ixgmin2=(ixomin2-1)*rft+1
3034 select case(direction_los)
3036 ixpmin1=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3037 ixpmax1=
node(pig2_,igrid)*rft*block_nx2
3038 ixpmin2=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3039 ixpmax2=
node(pig3_,igrid)*rft*block_nx3
3041 ixpmin1=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3042 ixpmax1=
node(pig3_,igrid)*rft*block_nx3
3043 ixpmin2=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3044 ixpmax2=
node(pig1_,igrid)*rft*block_nx1
3046 ixpmin1=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3047 ixpmax1=
node(pig1_,igrid)*rft*block_nx1
3048 ixpmin2=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3049 ixpmax2=
node(pig2_,igrid)*rft*block_nx2
3051 xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
3052 xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
3053 dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
3054 dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
3055 sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
3056 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
3058 deallocate(flux,dxb1,dxb2,dxb3,sxrg,xg1,xg2,dxg1,dxg2)
3065 integer,
intent(in) :: igrid,nXIF1,nXIF2
3066 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3067 double precision,
intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
3069 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
3071 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
3072 double precision :: xb^L,xd^D
3073 double precision,
allocatable :: flux(:^D&),v(:^D&),rho(:^D&),opacity(:^D&)
3074 double precision,
allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
3075 double precision,
allocatable :: EUVg(:,:),Fvg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
3076 integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
3077 double precision :: EUVt,Fvt,xc^L,xg^L,r2
3078 integer :: ixP^L,ixP^D
3079 integer :: direction_LOS
3089 ^d&ixomin^d=ixmlo^d\
3090 ^d&ixomax^d=ixmhi^d\
3091 ^d&iximin^d=
ixglo^d\
3092 ^d&iximax^d=
ixghi^d\
3096 allocate(flux(ixi^s),v(ixi^s),rho(ixi^s),opacity(ixi^s))
3097 allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
3098 dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
3099 dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
3100 dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
3113 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
3114 v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
3119 levelg=ps(igrid)%level
3123 select case(direction_los)
3134 allocate(euvg(nxg1,nxg2),fvg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
3141 select case(direction_los)
3143 do ix2=ixomin2,ixomax2
3144 ixgmin1=(ix2-1)*rft+1
3146 do ix3=ixomin3,ixomax3
3147 ixgmin2=(ix3-1)*rft+1
3151 do ix1=ixomin1,ixomax1
3155 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
3156 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
3160 do ix3=ixomin3,ixomax3
3161 ixgmin1=(ix3-1)*rft+1
3163 do ix1=ixomin1,ixomax1
3164 ixgmin2=(ix1-1)*rft+1
3168 do ix2=ixomin2,ixomax2
3172 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
3173 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
3177 do ix1=ixomin1,ixomax1
3178 ixgmin1=(ix1-1)*rft+1
3180 do ix2=ixomin2,ixomax2
3181 ixgmin2=(ix2-1)*rft+1
3185 do ix3=ixomin3,ixomax3
3189 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
3190 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
3201 select case(direction_los)
3203 ixgmin1=(ixomin2-1)*rft+1
3205 ixgmin2=(ixomin3-1)*rft+1
3208 ixgmin1=(ixomin3-1)*rft+1
3210 ixgmin2=(ixomin1-1)*rft+1
3213 ixgmin1=(ixomin1-1)*rft+1
3215 ixgmin2=(ixomin2-1)*rft+1
3219 select case(direction_los)
3221 ixpmin1=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3222 ixpmax1=
node(pig2_,igrid)*rft*block_nx2
3223 ixpmin2=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3224 ixpmax2=
node(pig3_,igrid)*rft*block_nx3
3226 ixpmin1=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3227 ixpmax1=
node(pig3_,igrid)*rft*block_nx3
3228 ixpmin2=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3229 ixpmax2=
node(pig1_,igrid)*rft*block_nx1
3231 ixpmin1=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3232 ixpmax1=
node(pig1_,igrid)*rft*block_nx1
3233 ixpmin2=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3234 ixpmax2=
node(pig2_,igrid)*rft*block_nx2
3236 xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
3237 xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
3238 dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
3239 dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
3240 euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
3241 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
3242 dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
3243 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
3245 deallocate(flux,v,opacity,dxb1,dxb2,dxb3,euvg,fvg,xg1,xg2,dxg1,dxg2)
3254 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),box_min(1:3),box_max(1:3)
3255 logical,
intent(out) :: hit
3256 double precision,
intent(out) :: t_enter,t_exit
3259 double precision :: t1,t2,td
3265 if (abs(ray_dir(idir))<=smalldouble)
then
3266 if (ray_origin(idir)<box_min(idir) .or. ray_origin(idir)>box_max(idir))
then
3271 t1=(box_min(idir)-ray_origin(idir))/ray_dir(idir)
3272 t2=(box_max(idir)-ray_origin(idir))/ray_dir(idir)
3278 t_enter=max(t_enter,t1)
3279 t_exit=min(t_exit,t2)
3280 if (t_enter>=t_exit)
then
3289 integer,
intent(in) :: ixI^L, ixO^L
3290 double precision,
intent(in) :: x(ixI^S,1:ndim),dx(ixI^S,1:ndim)
3291 double precision,
allocatable,
intent(out) :: xface1(:),xface2(:),xface3(:)
3295 allocate(xface1(ixomin1:ixomax1+1),xface2(ixomin2:ixomax2+1),xface3(ixomin3:ixomax3+1))
3299 do ix1=ixomin1,ixomax1
3300 xface1(ix1)=x(ix^d,1)-half*dx(ix^d,1)
3302 xface1(ixomax1+1)=x(ixomax1,ixomin2,ixomin3,1)+half*dx(ixomax1,ixomin2,ixomin3,1)
3306 do ix2=ixomin2,ixomax2
3307 xface2(ix2)=x(ix^d,2)-half*dx(ix^d,2)
3309 xface2(ixomax2+1)=x(ixomin1,ixomax2,ixomin3,2)+half*dx(ixomin1,ixomax2,ixomin3,2)
3313 do ix3=ixomin3,ixomax3
3314 xface3(ix3)=x(ix^d,3)-half*dx(ix^d,3)
3316 xface3(ixomax3+1)=x(ixomin1,ixomin2,ixomax3,3)+half*dx(ixomin1,ixomin2,ixomax3,3)
3320 integer,
intent(in) :: imin,imax
3321 double precision,
intent(in) :: pos,faces(imin:imax+1)
3323 integer :: ilo,ihi,imid
3325 if (pos<=faces(imin))
then
3329 if (pos>=faces(imax+1))
then
3336 do while (ihi-ilo>1)
3338 if (pos>=faces(imid))
then
3344 idx=min(imax,max(imin,ilo))
3348 integer,
intent(in) :: imin,imax,idx
3349 double precision,
intent(in) :: ray_origin_axis,ray_dir_axis,faces(imin:imax+1)
3350 integer,
intent(out) :: step
3351 double precision,
intent(out) :: tMax
3353 if (ray_dir_axis>zero)
then
3355 tmax=(faces(idx+1)-ray_origin_axis)/ray_dir_axis
3356 else if (ray_dir_axis<zero)
then
3358 tmax=(faces(idx)-ray_origin_axis)/ray_dir_axis
3366 integer,
intent(in) :: imin,imax,step
3367 double precision,
intent(in) :: ray_origin_axis,ray_dir_axis,faces(imin:imax+1)
3368 integer,
intent(inout) :: idx
3369 double precision,
intent(inout) :: tMax
3370 logical,
intent(out) :: done
3374 if (idx<imin .or. idx>imax)
then
3379 tmax=(faces(idx+1)-ray_origin_axis)/ray_dir_axis
3380 else if (step<0)
then
3381 tmax=(faces(idx)-ray_origin_axis)/ray_dir_axis
3388 ray_origin,xface1,xface2,xface3,t_enter,t_exit,EUVp,Dplp)
3389 integer,
intent(in) :: ixI^L, ixO^L
3390 double precision,
intent(in) :: source(ixI^S),sourcev(ixI^S)
3391 double precision,
intent(in) :: ray_origin(1:3)
3392 double precision,
intent(in) :: xface1(ixOmin1:ixOmax1+1),xface2(ixOmin2:ixOmax2+1),&
3393 xface3(ixOmin3:ixOmax3+1)
3394 double precision,
intent(in) :: t_enter,t_exit
3395 double precision,
intent(inout) :: EUVp,Dplp
3397 integer :: ix^D,step(1:3)
3398 double precision :: pos(1:3),tMax(1:3),tNow,tNext,ds_cm,epsRay
3401 if (t_exit<=t_enter)
return
3402 epsray=max(1.d-12,1.d-10*abs(t_exit-t_enter))
3403 pos=ray_origin+(t_enter+epsray)*
vec_los
3414 tnext=min(t_exit,tmax(1),tmax(2),tmax(3))
3415 if (tnext>tnow)
then
3416 ds_cm=(tnext-tnow)*unit_length
3417 if (si_unit) ds_cm=ds_cm*1.d2
3418 euvp=euvp+source(ix^d)*ds_cm
3419 dplp=dplp+sourcev(ix^d)*ds_cm
3422 if (tnow>=t_exit-epsray)
exit
3424 if (tmax(1)<=tnow+epsray)
then
3428 if (tmax(2)<=tnow+epsray)
then
3432 if (tmax(3)<=tnow+epsray)
then
3440 double precision,
allocatable,
intent(inout) :: segments(:,:)
3441 integer,
intent(inout) :: nseg,capacity
3442 integer,
intent(in) :: pixel_id
3443 double precision,
intent(in) :: tseg,jds,kds,jvds
3445 double precision,
allocatable :: tmp(:,:)
3446 integer :: new_capacity
3448 if (capacity<=0)
then
3450 allocate(segments(5,capacity))
3451 else if (nseg>=capacity)
then
3452 new_capacity=2*capacity
3453 allocate(tmp(5,new_capacity))
3454 tmp(:,1:capacity)=segments(:,1:capacity)
3455 call move_alloc(tmp,segments)
3456 capacity=new_capacity
3460 segments(1,nseg)=dble(pixel_id)
3461 segments(2,nseg)=tseg
3462 segments(3,nseg)=jds
3463 segments(4,nseg)=kds
3464 segments(5,nseg)=jvds
3468 pixel_id,ray_origin,xface1,xface2,xface3,t_enter,t_exit,&
3469 segments,nseg,capacity)
3470 integer,
intent(in) :: ixI^L, ixO^L
3471 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S),sourcev(ixI^S)
3472 integer,
intent(in) :: pixel_id
3473 double precision,
intent(in) :: ray_origin(1:3)
3474 double precision,
intent(in) :: xface1(ixOmin1:ixOmax1+1),xface2(ixOmin2:ixOmax2+1),&
3475 xface3(ixOmin3:ixOmax3+1)
3476 double precision,
intent(in) :: t_enter,t_exit
3477 double precision,
allocatable,
intent(inout) :: segments(:,:)
3478 integer,
intent(inout) :: nseg,capacity
3480 integer :: ix^D,step(1:3)
3481 double precision :: pos(1:3),tMax(1:3),tNow,tNext,ds_cm,epsRay,tseg
3482 double precision :: jds,kds,jvds
3485 if (t_exit<=t_enter)
return
3486 epsray=max(1.d-12,1.d-10*abs(t_exit-t_enter))
3487 pos=ray_origin+(t_enter+epsray)*
vec_los
3498 tnext=min(t_exit,tmax(1),tmax(2),tmax(3))
3499 if (tnext>tnow)
then
3500 ds_cm=(tnext-tnow)*unit_length
3501 if (si_unit) ds_cm=ds_cm*1.d2
3502 jds=source(ix^d)*ds_cm
3503 kds=opacity(ix^d)*ds_cm
3504 jvds=sourcev(ix^d)*ds_cm
3505 if (jds/=zero .or. kds/=zero .or. jvds/=zero)
then
3506 tseg=half*(tnow+tnext)
3511 if (tnow>=t_exit-epsray)
exit
3513 if (tmax(1)<=tnow+epsray)
then
3517 if (tmax(2)<=tnow+epsray)
then
3521 if (tmax(3)<=tnow+epsray)
then
3529 double precision,
intent(in) :: segments(:,:)
3530 integer,
intent(inout) :: idx(:)
3531 integer,
intent(in) :: nidx
3542 double precision,
intent(in) :: segments(:,:)
3543 integer,
intent(inout) :: idx(:)
3544 integer,
intent(in) :: ilo,ihi
3551 do while (j>=ilo .and. segments(2,idx(j))>segments(2,key))
3560 double precision,
intent(in) :: segments(:,:)
3561 integer,
intent(inout) :: idx(:)
3562 integer,
intent(in) :: ilo,ihi
3565 double precision :: pivot
3567 if (ihi-ilo<=32)
then
3574 pivot=segments(2,idx((ilo+ihi)/2))
3576 do while (segments(2,idx(i))<pivot)
3579 do while (segments(2,idx(j))>pivot)
3597 integer,
intent(in) :: pixel_id
3599 owner=mod(pixel_id-1,npe)
3603 double precision,
intent(in) :: segments(:,:)
3604 integer,
intent(in) :: is,nvars
3610 if (segments(iv,is)/=segments(iv,is) .or. abs(segments(iv,is))>=1.d90)
then
3617 subroutine cart_dda_block_pixel_range(box_min,box_max,nXIF1,nXIF2,xIF1,xIF2,ixPmin1,ixPmax1,ixPmin2,ixPmax2,has_pixels)
3618 double precision,
intent(in) :: box_min(1:3),box_max(1:3)
3619 integer,
intent(in) :: nXIF1,nXIF2
3620 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3621 integer,
intent(out) :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
3622 logical,
intent(out) :: has_pixels
3625 double precision :: vec_cor(1:3),xI_cor(1:2)
3626 double precision :: xmin1,xmax1,xmin2,xmax2,dx1,dx2
3629 if (i1==1) vec_cor(1)=box_min(1)
3630 if (i1==2) vec_cor(1)=box_max(1)
3632 if (i2==1) vec_cor(2)=box_min(2)
3633 if (i2==2) vec_cor(2)=box_max(2)
3635 if (i3==1) vec_cor(3)=box_min(3)
3636 if (i3==2) vec_cor(3)=box_max(3)
3638 if (i1==1 .and. i2==1 .and. i3==1)
then
3644 xmin1=min(xmin1,xi_cor(1))
3645 xmax1=max(xmax1,xi_cor(1))
3646 xmin2=min(xmin2,xi_cor(2))
3647 xmax2=max(xmax2,xi_cor(2))
3654 dx1=abs(xif1(2)-xif1(1))
3656 dx1=max(one,abs(xmax1-xmin1))
3659 dx2=abs(xif2(2)-xif2(1))
3661 dx2=max(one,abs(xmax2-xmin2))
3664 ixpmin1=max(1,floor((xmin1-xif1(1))/dx1)+1-1)
3665 ixpmax1=min(nxif1,ceiling((xmax1-xif1(1))/dx1)+1+1)
3666 ixpmin2=max(1,floor((xmin2-xif2(1))/dx2)+1-1)
3667 ixpmax2=min(nxif2,ceiling((xmax2-xif2(1))/dx2)+1+1)
3668 has_pixels=ixpmin1<=ixpmax1 .and. ixpmin2<=ixpmax2
3674 integer,
intent(in) :: nXIF1,nXIF2
3675 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3677 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
3679 integer :: ixO^L,ixO^D,ixI^L,ix^D
3680 integer :: iigrid,igrid,ixP1,ixP2,numSI,ixPmin1,ixPmax1,ixPmin2,ixPmax2
3681 double precision :: box_min(1:3),box_max(1:3),ray_origin(1:3)
3682 double precision :: t_enter,t_exit,vlos
3683 double precision :: profile_local(2),profile_global(2)
3684 logical :: hit,has_pixels
3685 double precision,
allocatable :: source(:^D&),sourcev(:^D&),rho(:^D&),opacity(:^D&)
3686 double precision,
allocatable :: xface1(:),xface2(:),xface3(:)
3687 double precision,
allocatable :: EUVs(:,:),Dpls(:,:)
3689 allocate(euvs(nxif1,nxif2),dpls(nxif1,nxif2))
3694 do iigrid=1,igridstail; igrid=igrids(iigrid);
3695 ^d&ixomin^d=ixmlo^d\
3696 ^d&ixomax^d=ixmhi^d\
3697 ^d&iximin^d=
ixglo^d\
3698 ^d&iximax^d=
ixghi^d\
3700 box_min(1)=
rnode(rpxmin1_,igrid)
3701 box_min(2)=
rnode(rpxmin2_,igrid)
3702 box_min(3)=
rnode(rpxmin3_,igrid)
3703 box_max(1)=
rnode(rpxmax1_,igrid)
3704 box_max(2)=
rnode(rpxmax2_,igrid)
3705 box_max(3)=
rnode(rpxmax3_,igrid)
3708 ixpmin1,ixpmax1,ixpmin2,ixpmax2,has_pixels)
3709 if (.not. has_pixels)
then
3710 deallocate(xface1,xface2,xface3)
3714 allocate(source(ixi^s),sourcev(ixi^s),rho(ixi^s),opacity(ixi^s))
3724 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
3725 do ix1=ixomin1,ixomax1
3726 do ix2=ixomin2,ixomax2
3727 do ix3=ixomin3,ixomax3
3728 if (rho(ix^d)>smalldouble)
then
3729 vlos=(ps(igrid)%w(ix^d,iw_mom(1))*
vec_los(1)+&
3730 ps(igrid)%w(ix^d,iw_mom(2))*
vec_los(2)+&
3731 ps(igrid)%w(ix^d,iw_mom(3))*
vec_los(3))/rho(ix^d)
3732 sourcev(ix^d)=source(ix^d)*vlos
3738 deallocate(rho,opacity)
3740 do ixp1=ixpmin1,ixpmax1
3741 do ixp2=ixpmin2,ixpmax2
3743 profile_local(1)=profile_local(1)+one
3746 profile_local(2)=profile_local(2)+one
3748 ray_origin,xface1,xface2,xface3,t_enter,t_exit,euvs(ixp1,ixp2),dpls(ixp1,ixp2))
3753 deallocate(source,sourcev,xface1,xface2,xface3)
3757 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
3758 call mpi_allreduce(dpls,dpl,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
3759 call mpi_allreduce(profile_local,profile_global,2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
3761 write(*,
'(a,2(es12.5,1x))')
' cart_dda thin profile ray_tests ray_hits: ',profile_global
3763 deallocate(euvs,dpls)
3769 integer,
intent(in) :: nXIF1,nXIF2
3770 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3772 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
3773 double precision,
intent(out) :: Tau(nXIF1,nXIF2),EUVthin(nXIF1,nXIF2)
3775 integer,
parameter :: nSegVars=5
3776 integer :: ixO^L,ixO^D,ixI^L,ix^D
3777 integer :: iigrid,igrid,ixP1,ixP2,ipix,ipixStart,ipixEnd,nPixBatch,pixel_id
3778 integer :: nseg,capacity,totalCount,totalSeg,ipe,is,iseg,nidx,owner,isegDest,nsegBefore
3779 integer :: ixGlobal,iyGlobal,ixPmin1,ixPmax1,ixPmin2,ixPmax2,iFirst,iLast,iLocal
3780 integer :: nPixBatchTarget
3781 integer :: maxSegBatchTarget,maxSegCommTarget,maxNsegBatch,nPixTotal
3782 integer :: maxOwnerSegCount,maxOwnerSegCountLocal,segOffset,recvFill,totalRoundCount,totalRoundSeg
3783 integer,
allocatable :: sendCounts(:),recvCounts(:),sendDispls(:),recvDispls(:)
3784 integer,
allocatable :: roundSendCounts(:),roundRecvCounts(:)
3785 integer,
allocatable :: roundSendDispls(:),roundRecvDispls(:)
3786 integer,
allocatable :: ownerSegCounts(:),ownerOffsets(:),idx(:)
3787 integer,
allocatable :: bucketCounts(:),bucketOffsets(:),bucketFill(:)
3788 double precision :: ray_origin(1:3)
3789 double precision :: t_enter,t_exit,vlos,atten
3790 double precision :: profile_local(5),profile_global(5),profile_batch(5)
3791 logical :: hit,has_pixels,batchAccepted,batchReduced
3792 double precision,
allocatable :: rho(:^D&)
3793 double precision,
allocatable :: segments(:,:),segments_send(:,:),segments_recv(:,:)
3794 double precision,
allocatable :: segments_recv_round(:,:)
3795 double precision,
allocatable :: image_reduce(:,:)
3803 allocate(sendcounts(0:
npe-1),recvcounts(0:
npe-1),senddispls(0:
npe-1),recvdispls(0:
npe-1))
3804 allocate(roundsendcounts(0:
npe-1),roundrecvcounts(0:
npe-1))
3805 allocate(roundsenddispls(0:
npe-1),roundrecvdispls(0:
npe-1))
3806 allocate(ownersegcounts(0:
npe-1),owneroffsets(0:
npe-1))
3807 allocate(cache(igridstail))
3809 allocate(bucketcounts(npixbatchtarget),bucketoffsets(npixbatchtarget+1),&
3810 bucketfill(npixbatchtarget))
3812 do iigrid=1,igridstail; igrid=igrids(iigrid);
3813 ^d&ixomin^d=ixmlo^d\
3814 ^d&ixomax^d=ixmhi^d\
3815 ^d&iximin^d=
ixglo^d\
3816 ^d&iximax^d=
ixghi^d\
3818 cache(iigrid)%igrid=igrid
3819 allocate(cache(iigrid)%source(ixi^s),cache(iigrid)%opacity(ixi^s),&
3820 cache(iigrid)%sourcev(ixi^s),rho(ixi^s))
3821 cache(iigrid)%source=zero
3822 cache(iigrid)%opacity=zero
3823 cache(iigrid)%sourcev=zero
3826 cache(iigrid)%source,cache(iigrid)%opacity)
3828 call get_euv(
wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,cache(iigrid)%source)
3831 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
3832 do ix1=ixomin1,ixomax1
3833 do ix2=ixomin2,ixomax2
3834 do ix3=ixomin3,ixomax3
3835 if (rho(ix^d)>smalldouble)
then
3836 vlos=(ps(igrid)%w(ix^d,iw_mom(1))*
vec_los(1)+&
3837 ps(igrid)%w(ix^d,iw_mom(2))*
vec_los(2)+&
3838 ps(igrid)%w(ix^d,iw_mom(3))*
vec_los(3))/rho(ix^d)
3839 cache(iigrid)%sourcev(ix^d)=cache(iigrid)%source(ix^d)*vlos
3846 cache(iigrid)%box_min(1)=
rnode(rpxmin1_,igrid)
3847 cache(iigrid)%box_min(2)=
rnode(rpxmin2_,igrid)
3848 cache(iigrid)%box_min(3)=
rnode(rpxmin3_,igrid)
3849 cache(iigrid)%box_max(1)=
rnode(rpxmax1_,igrid)
3850 cache(iigrid)%box_max(2)=
rnode(rpxmax2_,igrid)
3851 cache(iigrid)%box_max(3)=
rnode(rpxmax3_,igrid)
3853 cache(iigrid)%xface1,cache(iigrid)%xface2,&
3854 cache(iigrid)%xface3)
3856 nxif1,nxif2,xif1,xif2,cache(iigrid)%ixPmin1,cache(iigrid)%ixPmax1,&
3857 cache(iigrid)%ixPmin2,cache(iigrid)%ixPmax2,cache(iigrid)%has_pixels)
3860 npixtotal=nxif1*nxif2
3862 do while (ipixstart<=npixtotal)
3863 ipixend=min(nxif1*nxif2,ipixstart+npixbatchtarget-1)
3864 npixbatch=ipixend-ipixstart+1
3865 batchaccepted=.false.
3866 batchreduced=.false.
3868 do while (.not. batchaccepted)
3873 do iigrid=1,igridstail; igrid=igrids(iigrid);
3874 ^d&ixomin^d=ixmlo^d\
3875 ^d&ixomax^d=ixmhi^d\
3876 ^d&iximin^d=
ixglo^d\
3877 ^d&iximax^d=
ixghi^d\
3879 ixpmin1=cache(iigrid)%ixPmin1
3880 ixpmax1=cache(iigrid)%ixPmax1
3881 ixpmin2=cache(iigrid)%ixPmin2
3882 ixpmax2=cache(iigrid)%ixPmax2
3883 has_pixels=cache(iigrid)%has_pixels
3884 if (.not. has_pixels) cycle
3886 do ixp2=ixpmin2,ixpmax2
3887 ifirst=max(ipixstart,(ixp2-1)*nxif1+ixpmin1)
3888 ilast=min(ipixend,(ixp2-1)*nxif1+ixpmax1)
3889 if (ifirst>ilast) cycle
3890 do ipix=ifirst,ilast
3891 ixp1=1+mod(ipix-1,nxif1)
3893 profile_batch(1)=profile_batch(1)+one
3895 cache(iigrid)%box_max,hit,t_enter,t_exit)
3897 profile_batch(2)=profile_batch(2)+one
3900 cache(iigrid)%opacity,cache(iigrid)%sourcev,&
3901 ipix,ray_origin,cache(iigrid)%xface1,&
3902 cache(iigrid)%xface2,cache(iigrid)%xface3,&
3904 segments,nseg,capacity)
3905 profile_batch(3)=profile_batch(3)+dble(nseg-nsegbefore)
3911 call mpi_allreduce(nseg,maxnsegbatch,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
3912 if (maxnsegbatch>maxsegbatchtarget .and. npixbatch>1)
then
3913 npixbatch=max(1,npixbatch/2)
3914 ipixend=ipixstart+npixbatch-1
3915 if (
allocated(segments))
deallocate(segments)
3918 batchaccepted=.true.
3922 profile_local=profile_local+profile_batch
3924 write(*,
'(a,3(i0,1x))')
' cart_dda thick adaptive batch: ',&
3925 ipixstart,ipixend,maxnsegbatch
3928 if (.not.
allocated(segments))
then
3930 allocate(segments(nsegvars,capacity))
3935 ownersegcounts(owner)=ownersegcounts(owner)+1
3937 sendcounts=nsegvars*ownersegcounts
3940 senddispls(ipe)=senddispls(ipe-1)+sendcounts(ipe-1)
3943 allocate(segments_send(nsegvars,max(1,nseg)))
3947 isegdest=senddispls(owner)/nsegvars+owneroffsets(owner)+1
3948 segments_send(:,isegdest)=segments(:,is)
3949 owneroffsets(owner)=owneroffsets(owner)+1
3952 call mpi_alltoall(sendcounts,1,mpi_integer,recvcounts,1,mpi_integer,
icomm,
ierrmpi)
3955 recvdispls(ipe)=recvdispls(ipe-1)+recvcounts(ipe-1)
3957 totalcount=sum(recvcounts)
3958 totalseg=totalcount/nsegvars
3959 profile_local(4)=profile_local(4)+dble(totalcount)
3960 allocate(segments_recv(nsegvars,max(1,totalseg)))
3963 maxownersegcountlocal=maxval(ownersegcounts)
3964 call mpi_allreduce(maxownersegcountlocal,maxownersegcount,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
3965 do segoffset=0,maxownersegcount-1,maxsegcommtarget
3967 roundsenddispls=senddispls
3969 if (ownersegcounts(ipe)>segoffset)
then
3970 roundsendcounts(ipe)=nsegvars*min(maxsegcommtarget,ownersegcounts(ipe)-segoffset)
3971 roundsenddispls(ipe)=senddispls(ipe)+nsegvars*segoffset
3975 call mpi_alltoall(roundsendcounts,1,mpi_integer,roundrecvcounts,1,mpi_integer,
icomm,
ierrmpi)
3976 roundrecvdispls(0)=0
3978 roundrecvdispls(ipe)=roundrecvdispls(ipe-1)+roundrecvcounts(ipe-1)
3980 totalroundcount=sum(roundrecvcounts)
3981 totalroundseg=totalroundcount/nsegvars
3982 allocate(segments_recv_round(nsegvars,max(1,totalroundseg)))
3984 call mpi_alltoallv(segments_send,roundsendcounts,roundsenddispls,mpi_double_precision,&
3985 segments_recv_round,roundrecvcounts,roundrecvdispls,&
3988 if (totalroundseg>0)
then
3989 segments_recv(:,recvfill+1:recvfill+totalroundseg)=segments_recv_round(:,1:totalroundseg)
3990 recvfill=recvfill+totalroundseg
3992 deallocate(segments_recv_round)
3995 if (recvfill/=totalseg)
call mpistop(
"cart_dda thick segmented receive mismatch")
3997 if (totalseg>0)
then
3998 allocate(idx(totalseg))
3999 bucketcounts(1:npixbatch)=0
4002 ipix=nint(segments_recv(1,is))
4004 ilocal=ipix-ipixstart+1
4005 bucketcounts(ilocal)=bucketcounts(ilocal)+1
4011 do ilocal=1,npixbatch
4012 bucketoffsets(ilocal+1)=bucketoffsets(ilocal)+bucketcounts(ilocal)
4014 bucketfill(1:npixbatch)=bucketoffsets(1:npixbatch)
4017 ipix=nint(segments_recv(1,is))
4019 ilocal=ipix-ipixstart+1
4020 idx(bucketfill(ilocal))=is
4021 bucketfill(ilocal)=bucketfill(ilocal)+1
4026 do ipix=ipixstart,ipixend
4028 ilocal=ipix-ipixstart+1
4029 nidx=bucketcounts(ilocal)
4031 profile_local(5)=profile_local(5)+dble(nidx)*dble(nidx)
4033 ixglobal=1+mod(ipix-1,nxif1)
4034 iyglobal=1+(ipix-1)/nxif1
4035 do iseg=bucketoffsets(ilocal),bucketoffsets(ilocal+1)-1
4037 euvthin(ixglobal,iyglobal)=euvthin(ixglobal,iyglobal)+segments_recv(3,is)
4039 euv(ixglobal,iyglobal)=euv(ixglobal,iyglobal)+atten*segments_recv(3,is)
4040 dpl(ixglobal,iyglobal)=dpl(ixglobal,iyglobal)+atten*segments_recv(5,is)
4041 tau(ixglobal,iyglobal)=tau(ixglobal,iyglobal)+max(zero,segments_recv(4,is))
4048 deallocate(segments_send,segments_recv)
4049 if (
allocated(segments))
deallocate(segments)
4053 do iigrid=1,igridstail
4054 if (
allocated(cache(iigrid)%source))
deallocate(cache(iigrid)%source)
4055 if (
allocated(cache(iigrid)%opacity))
deallocate(cache(iigrid)%opacity)
4056 if (
allocated(cache(iigrid)%sourcev))
deallocate(cache(iigrid)%sourcev)
4057 if (
allocated(cache(iigrid)%xface1))
deallocate(cache(iigrid)%xface1)
4058 if (
allocated(cache(iigrid)%xface2))
deallocate(cache(iigrid)%xface2)
4059 if (
allocated(cache(iigrid)%xface3))
deallocate(cache(iigrid)%xface3)
4062 deallocate(sendcounts,recvcounts,senddispls,recvdispls,roundsendcounts,roundrecvcounts,&
4063 roundsenddispls,roundrecvdispls,ownersegcounts,owneroffsets,bucketcounts,&
4064 bucketoffsets,bucketfill)
4065 allocate(image_reduce(nxif1,nxif2))
4066 call mpi_allreduce(euv,image_reduce,nxif1*nxif2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4068 call mpi_allreduce(dpl,image_reduce,nxif1*nxif2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4070 call mpi_allreduce(tau,image_reduce,nxif1*nxif2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4072 call mpi_allreduce(euvthin,image_reduce,nxif1*nxif2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4073 euvthin=image_reduce
4074 deallocate(image_reduce)
4075 call mpi_allreduce(profile_local,profile_global,5,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4077 write(*,
'(a,5(es12.5,1x))') &
4078 ' cart_dda thick profile: ',profile_global
4085 integer,
intent(in) :: nXIF1,nXIF2
4087 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
4088 double precision,
intent(out) :: Tau(nXIF1,nXIF2),EUVthin(nXIF1,nXIF2)
4090 integer :: ixO^L,ixO^D,ixI^L,ix^D
4091 integer :: iigrid,igrid,levelg,rft,direction_LOS,nLOS,numSeg,nLayerVars,nLayerSeg
4092 integer :: ixP1,ixP2,ixL,iSub1,iSub2,relL
4093 integer :: nLosBatch,nBatch,iBatch,ixLstart,ixLend,ixLgridStart,ixLgridEnd
4094 integer :: ixPmin1,ixPmin2
4095 double precision :: ds_cm,jds,kds,jvds,atten,layerBytes,targetBytes
4096 double precision,
allocatable :: rho(:^D&)
4097 double precision,
allocatable :: layer_ds(:,:,:,:),layer_all(:,:,:,:)
4112 if (nxif1>huge(numseg)/max(1,nxif2) .or. nxif1*nxif2>huge(numseg)/nlayervars)
then
4113 call mpistop(
"thick EUV layer buffer is too large for one MPI reduction")
4115 nlayerseg=nxif1*nxif2*nlayervars
4116 targetbytes=256.d0*1024.d0*1024.d0
4117 layerbytes=dble(nlayerseg)*8.d0*2.d0
4118 nlosbatch=max(1,min(16,int(targetbytes/max(one,layerbytes))))
4119 if (nlayerseg>huge(numseg)/nlosbatch)
then
4120 call mpistop(
"thick EUV batched layer buffer is too large for one MPI reduction")
4123 allocate(cache(igridstail))
4124 do iigrid=1,igridstail; igrid=igrids(iigrid);
4125 ^d&ixomin^d=ixmlo^d\
4126 ^d&ixomax^d=ixmhi^d\
4127 ^d&iximin^d=
ixglo^d\
4128 ^d&iximax^d=
ixghi^d\
4130 cache(iigrid)%igrid=igrid
4131 levelg=ps(igrid)%level
4133 cache(iigrid)%level=levelg
4134 cache(iigrid)%rft=rft
4136 select case(direction_los)
4138 cache(iigrid)%los_min=(
node(pig1_,igrid)-1)*rft*block_nx1+1
4139 cache(iigrid)%los_max=
node(pig1_,igrid)*rft*block_nx1
4141 cache(iigrid)%los_min=(
node(pig2_,igrid)-1)*rft*block_nx2+1
4142 cache(iigrid)%los_max=
node(pig2_,igrid)*rft*block_nx2
4144 cache(iigrid)%los_min=(
node(pig3_,igrid)-1)*rft*block_nx3+1
4145 cache(iigrid)%los_max=
node(pig3_,igrid)*rft*block_nx3
4148 allocate(cache(iigrid)%source(ixi^s),cache(iigrid)%opacity(ixi^s),&
4149 cache(iigrid)%sourcev(ixi^s),rho(ixi^s))
4150 cache(iigrid)%source=zero
4151 cache(iigrid)%opacity=zero
4152 cache(iigrid)%sourcev=zero
4155 cache(iigrid)%source,cache(iigrid)%opacity)
4157 call get_euv(
wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,cache(iigrid)%source)
4160 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
4161 cache(iigrid)%sourcev(ixo^s)=cache(iigrid)%source(ixo^s)*&
4162 (-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s))
4168 allocate(layer_ds(nxif1,nxif2,nlayervars,nlosbatch),layer_all(nxif1,nxif2,nlayervars,nlosbatch))
4174 do ixlstart=1,nlos,nlosbatch
4175 ixlend=min(nlos,ixlstart+nlosbatch-1)
4176 nbatch=ixlend-ixlstart+1
4177 layer_ds(:,:,:,1:nbatch)=zero
4179 do iigrid=1,igridstail
4180 ixlgridstart=max(ixlstart,cache(iigrid)%los_min)
4181 ixlgridend=min(ixlend,cache(iigrid)%los_max)
4182 if (ixlgridstart>ixlgridend) cycle
4183 igrid=cache(iigrid)%igrid
4184 rft=cache(iigrid)%rft
4185 ^d&ixomin^d=ixmlo^d\
4186 ^d&ixomax^d=ixmhi^d\
4187 ^d&iximin^d=
ixglo^d\
4188 ^d&iximax^d=
ixghi^d\
4190 do ixl=ixlgridstart,ixlgridend
4191 ibatch=ixl-ixlstart+1
4192 rell=ixl-cache(iigrid)%los_min
4194 select case(direction_los)
4196 ix1=ixomin1+rell/rft
4197 do ix2=ixomin2,ixomax2
4198 ixpmin1=(
node(pig2_,igrid)-1)*rft*block_nx2+(ix2-ixomin2)*rft+1
4199 do ix3=ixomin3,ixomax3
4200 ixpmin2=(
node(pig3_,igrid)-1)*rft*block_nx3+(ix3-ixomin3)*rft+1
4203 jds=cache(iigrid)%source(ix^d)*ds_cm
4204 kds=cache(iigrid)%opacity(ix^d)*ds_cm
4205 jvds=cache(iigrid)%sourcev(ix^d)*ds_cm
4210 layer_ds(ixp1,ixp2,1,ibatch)=layer_ds(ixp1,ixp2,1,ibatch)+jds
4211 layer_ds(ixp1,ixp2,2,ibatch)=layer_ds(ixp1,ixp2,2,ibatch)+kds
4212 layer_ds(ixp1,ixp2,3,ibatch)=layer_ds(ixp1,ixp2,3,ibatch)+jvds
4218 ix2=ixomin2+rell/rft
4219 do ix3=ixomin3,ixomax3
4220 ixpmin1=(
node(pig3_,igrid)-1)*rft*block_nx3+(ix3-ixomin3)*rft+1
4221 do ix1=ixomin1,ixomax1
4222 ixpmin2=(
node(pig1_,igrid)-1)*rft*block_nx1+(ix1-ixomin1)*rft+1
4225 jds=cache(iigrid)%source(ix^d)*ds_cm
4226 kds=cache(iigrid)%opacity(ix^d)*ds_cm
4227 jvds=cache(iigrid)%sourcev(ix^d)*ds_cm
4232 layer_ds(ixp1,ixp2,1,ibatch)=layer_ds(ixp1,ixp2,1,ibatch)+jds
4233 layer_ds(ixp1,ixp2,2,ibatch)=layer_ds(ixp1,ixp2,2,ibatch)+kds
4234 layer_ds(ixp1,ixp2,3,ibatch)=layer_ds(ixp1,ixp2,3,ibatch)+jvds
4240 ix3=ixomin3+rell/rft
4241 do ix1=ixomin1,ixomax1
4242 ixpmin1=(
node(pig1_,igrid)-1)*rft*block_nx1+(ix1-ixomin1)*rft+1
4243 do ix2=ixomin2,ixomax2
4244 ixpmin2=(
node(pig2_,igrid)-1)*rft*block_nx2+(ix2-ixomin2)*rft+1
4247 jds=cache(iigrid)%source(ix^d)*ds_cm
4248 kds=cache(iigrid)%opacity(ix^d)*ds_cm
4249 jvds=cache(iigrid)%sourcev(ix^d)*ds_cm
4254 layer_ds(ixp1,ixp2,1,ibatch)=layer_ds(ixp1,ixp2,1,ibatch)+jds
4255 layer_ds(ixp1,ixp2,2,ibatch)=layer_ds(ixp1,ixp2,2,ibatch)+kds
4256 layer_ds(ixp1,ixp2,3,ibatch)=layer_ds(ixp1,ixp2,3,ibatch)+jvds
4265 numseg=nlayerseg*nbatch
4266 call mpi_allreduce(layer_ds,layer_all,numseg,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4272 euvthin(ixp1,ixp2)=euvthin(ixp1,ixp2)+layer_all(ixp1,ixp2,1,ibatch)
4274 euv(ixp1,ixp2)=euv(ixp1,ixp2)+atten*layer_all(ixp1,ixp2,1,ibatch)
4275 dpl(ixp1,ixp2)=dpl(ixp1,ixp2)+atten*layer_all(ixp1,ixp2,3,ibatch)
4276 tau(ixp1,ixp2)=tau(ixp1,ixp2)+layer_all(ixp1,ixp2,2,ibatch)
4283 do iigrid=1,igridstail
4284 if (
allocated(cache(iigrid)%source))
deallocate(cache(iigrid)%source)
4285 if (
allocated(cache(iigrid)%opacity))
deallocate(cache(iigrid)%opacity)
4286 if (
allocated(cache(iigrid)%sourcev))
deallocate(cache(iigrid)%sourcev)
4288 deallocate(cache,layer_ds,layer_all)
4297 integer,
intent(in) :: ixI^L, ixO^L
4298 double precision,
intent(in) :: x(ixI^S,1:ndim),dx(ixI^S,1:ndim)
4299 double precision,
allocatable,
intent(out) :: rface(:),thetaface(:),phiface(:)
4300 integer :: ix1,ix2,ix3
4302 allocate(rface(ixomin1:ixomax1+1),thetaface(ixomin2:ixomax2+1),phiface(ixomin3:ixomax3+1))
4303 do ix1=ixomin1,ixomax1
4304 rface(ix1)=x(ix1,ixomin2,ixomin3,1)-half*dx(ix1,ixomin2,ixomin3,1)
4306 rface(ixomax1+1)=x(ixomax1,ixomin2,ixomin3,1)+half*dx(ixomax1,ixomin2,ixomin3,1)
4307 do ix2=ixomin2,ixomax2
4308 thetaface(ix2)=x(ixomin1,ix2,ixomin3,2)-half*dx(ixomin1,ix2,ixomin3,2)
4310 thetaface(ixomax2+1)=x(ixomin1,ixomax2,ixomin3,2)+half*dx(ixomin1,ixomax2,ixomin3,2)
4311 do ix3=ixomin3,ixomax3
4312 phiface(ix3)=x(ixomin1,ixomin2,ix3,3)-half*dx(ixomin1,ixomin2,ix3,3)
4314 phiface(ixomax3+1)=x(ixomin1,ixomin2,ixomax3,3)+half*dx(ixomin1,ixomin2,ixomax3,3)
4318 double precision,
allocatable,
intent(inout) :: tvals(:)
4319 integer,
intent(inout) :: nt,capacity
4320 double precision,
intent(in) :: t
4322 double precision,
allocatable :: tmp(:)
4324 if (t /= t .or. abs(t)>1.d90)
return
4325 if (.not.
allocated(tvals))
then
4327 allocate(tvals(capacity))
4328 else if (nt>=capacity)
then
4329 allocate(tmp(capacity))
4332 allocate(tvals(2*capacity))
4333 tvals(1:capacity)=tmp
4342 double precision,
intent(inout) :: tvals(:)
4343 integer,
intent(inout) :: nt
4344 double precision,
intent(in) :: t
4346 if (t /= t .or. abs(t)>1.d90)
return
4347 if (nt>=
size(tvals))
return
4353 double precision,
intent(inout) :: tvals(:)
4354 integer,
intent(inout) :: nt
4357 double precision :: key,epsT
4363 do while (j>=1 .and. tvals(j)>key)
4369 epst=max(1.d-12,1.d-10*max(one,abs(tvals(nt)-tvals(1))))
4372 if (abs(tvals(i)-tvals(nout))>epst)
then
4374 tvals(nout)=tvals(i)
4381 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),rface
4382 double precision,
allocatable,
intent(inout) :: tvals(:)
4383 integer,
intent(inout) :: nt,capacity
4385 double precision :: aa,bb,cc,disc,root
4388 bb=2.d0*sum(ray_origin*ray_dir)
4389 cc=sum(ray_origin**2)-rface**2
4390 disc=bb**2-4.d0*aa*cc
4391 if (disc<zero)
return
4392 root=sqrt(max(zero,disc))
4393 call sph_add_t(tvals,nt,capacity,(-bb-root)/(2.d0*aa))
4394 call sph_add_t(tvals,nt,capacity,(-bb+root)/(2.d0*aa))
4398 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),thetaface
4399 double precision,
allocatable,
intent(inout) :: tvals(:)
4400 integer,
intent(inout) :: nt,capacity
4402 double precision :: cth,aa,bb,cc,disc,root
4405 aa=ray_dir(3)**2-cth**2*sum(ray_dir**2)
4406 bb=2.d0*(ray_origin(3)*ray_dir(3)-cth**2*sum(ray_origin*ray_dir))
4407 cc=ray_origin(3)**2-cth**2*sum(ray_origin**2)
4408 if (abs(aa)<1.d-14)
then
4409 if (abs(bb)>1.d-14)
call sph_add_t(tvals,nt,capacity,-cc/bb)
4412 disc=bb**2-4.d0*aa*cc
4413 if (disc<zero)
return
4414 root=sqrt(max(zero,disc))
4415 call sph_add_t(tvals,nt,capacity,(-bb-root)/(2.d0*aa))
4416 call sph_add_t(tvals,nt,capacity,(-bb+root)/(2.d0*aa))
4420 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),phiface
4421 double precision,
allocatable,
intent(inout) :: tvals(:)
4422 integer,
intent(inout) :: nt,capacity
4424 double precision :: normal(1:3),denom,numer
4426 normal(1)=-sin(phiface)
4427 normal(2)=cos(phiface)
4429 denom=sum(normal*ray_dir)
4430 if (abs(denom)<1.d-14)
return
4431 numer=sum(normal*ray_origin)
4432 call sph_add_t(tvals,nt,capacity,-numer/denom)
4436 double precision,
intent(in) :: pos(1:3)
4437 double precision,
intent(out) :: sph(1:3)
4439 sph(1)=sqrt(sum(pos**2))
4440 if (sph(1)>zero)
then
4441 sph(2)=acos(max(-one,min(one,pos(3)/sph(1))))
4445 sph(3)=atan2(pos(2),pos(1))
4449 integer,
intent(in) :: imin,imax
4450 double precision,
intent(in) ::
value,faces(imin:imax+1)
4452 integer :: ilo,ihi,imid
4455 if (
value<faces(imin)-1.d-12 .or.
value>faces(imax+1)+1.d-12)
return
4456 if (
value<=faces(imin))
then
4460 if (
value>=faces(imax+1))
then
4467 do while (ihi-ilo>1)
4469 if (
value>=faces(imid))
then
4475 idx=min(imax,max(imin,ilo))
4479 integer,
intent(in) :: imin,imax
4480 double precision,
intent(in) ::
value,faces(imin:imax+1)
4482 integer :: ilo,ihi,imid
4485 if (
value>faces(imin)+1.d-12 .or.
value<faces(imax+1)-1.d-12)
return
4486 if (
value>=faces(imin))
then
4490 if (
value<=faces(imax+1))
then
4497 do while (ihi-ilo>1)
4499 if (
value<=faces(imid))
then
4505 idx=min(imax,max(imin,ilo))
4509 double precision,
intent(in) :: pos(1:3)
4510 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4511 phiface(ixOmin3:ixOmax3+1)
4512 integer,
intent(in) :: ixO^L
4513 integer,
intent(out) :: ix1,ix2,ix3
4514 logical,
intent(out) :: inside
4516 double precision :: sph(1:3),phi
4520 if (phi<phiface(ixomin3)-1.d-12) phi=phi+2.d0*dpi
4521 if (phi>phiface(ixomax3+1)+1.d-12) phi=phi-2.d0*dpi
4525 inside=ix1>0 .and. ix2>0 .and. ix3>0
4530 double precision,
intent(in) :: pos(1:3),ximg1,ximg2
4532 double precision :: dotp,rc,rthick,rloc
4535 rc=sqrt(ximg1**2+ximg2**2)
4536 rloc=sqrt(sum(pos**2))
4539 if (dotp>=
zero)
then
4540 if (rc<=rthick) visible=.false.
4542 if (rloc<=rthick) visible=.false.
4547 ixPmin1,ixPmax1,ixPmin2,ixPmax2,has_pixels)
4548 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4549 phiface(ixOmin3:ixOmax3+1)
4550 integer,
intent(in) :: ixO^L,nXI1,nXI2
4551 double precision,
intent(in) :: xI1(nXI1),xI2(nXI2),dxI
4552 integer,
intent(out) :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
4553 logical,
intent(out) :: has_pixels
4555 integer,
parameter :: nsample=5
4557 double precision :: sph(1:3),xcent(1:2)
4558 double precision :: xmin1,xmax1,xmin2,xmax2
4559 double precision :: wr,wt,wp,pad
4567 wr=dble(ir)/dble(nsample-1)
4568 sph(1)=(one-wr)*rface(ixomin1)+wr*rface(ixomax1+1)
4570 wt=dble(it)/dble(nsample-1)
4571 sph(2)=(one-wt)*thetaface(ixomin2)+wt*thetaface(ixomax2+1)
4573 wp=dble(ip)/dble(nsample-1)
4574 if (ir/=0 .and. ir/=nsample-1 .and. it/=0 .and. it/=nsample-1 .and. &
4575 ip/=0 .and. ip/=nsample-1) cycle
4576 sph(3)=(one-wp)*phiface(ixomin3)+wp*phiface(ixomax3+1)
4578 xmin1=min(xmin1,xcent(1))
4579 xmax1=max(xmax1,xcent(1))
4580 xmin2=min(xmin2,xcent(2))
4581 xmax2=max(xmax2,xcent(2))
4590 ixpmin1=max(1,floor((xmin1-(xi1(1)-half*dxi))/dxi)+1)
4591 ixpmax1=min(nxi1,ceiling((xmax1-(xi1(1)-half*dxi))/dxi))
4592 ixpmin2=max(1,floor((xmin2-(xi2(1)-half*dxi))/dxi)+1)
4593 ixpmax2=min(nxi2,ceiling((xmax2-(xi2(1)-half*dxi))/dxi))
4594 has_pixels=ixpmin1<=ixpmax1 .and. ixpmin2<=ixpmax2
4598 rface,thetaface,phiface,EUVp)
4599 integer,
intent(in) :: ixI^L,ixO^L
4600 double precision,
intent(in) :: source(ixI^S),ray_origin(1:3),ximg1,ximg2
4601 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4602 phiface(ixOmin3:ixOmax3+1)
4603 double precision,
intent(inout) :: EUVp
4605 integer :: nt,capacity,i,ix^D
4606 double precision,
allocatable :: tvals(:)
4607 double precision :: posMid(1:3),ds_cm,tMid,t0,t1
4612 do ix1=ixomin1,ixomax1+1
4615 do ix2=ixomin2,ixomax2+1
4618 do ix3=ixomin3,ixomax3+1
4622 if (
allocated(tvals))
deallocate(tvals)
4632 posmid=ray_origin+tmid*
vec_los
4634 call sph_locate_cell(posmid,rface,thetaface,phiface,ixo^l,ix1,ix2,ix3,inside)
4635 if (.not. inside) cycle
4636 ds_cm=(t1-t0)*unit_length
4637 if (si_unit) ds_cm=ds_cm*1.d2
4638 euvp=euvp+source(ix^d)*ds_cm
4644 pixel_id,ray_origin,ximg1,ximg2,&
4645 rface,thetaface,phiface,rface2,&
4646 theta_cos,phi_sin,phi_cos,&
4647 segments,nseg,capacity)
4650 integer,
intent(in) :: ixI^L,ixO^L,pixel_id
4651 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S)
4652 double precision,
intent(in) :: ray_origin(1:3),ximg1,ximg2
4653 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4654 phiface(ixOmin3:ixOmax3+1)
4655 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4656 phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
4657 double precision,
allocatable,
intent(inout) :: segments(:,:)
4658 integer,
intent(inout) :: nseg,capacity
4660 integer :: nt,i,ix^D
4661 double precision :: tvals(2*(ixOmax1-ixOmin1+2)+2*(ixOmax2-ixOmin2+2)+&
4662 (ixOmax3-ixOmin3+2))
4663 double precision :: posMid(1:3),ds_cm,tMid,t0,t1,jds,kds
4664 double precision :: dir2,origin2,odotd,aa,bb,cc,disc,root,cth2
4665 double precision :: denom,numer,r2,mu,phi,dotp,rthick2,rc2
4670 origin2=sum(ray_origin**2)
4672 rthick2=(r_opt_thick*
const_rsun/unit_length)**2
4673 rc2=ximg1**2+ximg2**2
4675 do ix1=ixomin1,ixomax1+1
4678 cc=origin2-rface2(ix1)
4679 disc=bb**2-4.d0*aa*cc
4680 if (disc>=
zero)
then
4681 root=sqrt(max(
zero,disc))
4686 do ix2=ixomin2,ixomax2+1
4687 cth2=theta_cos(ix2)**2
4689 bb=2.d0*(ray_origin(3)*
vec_los(3)-cth2*odotd)
4690 cc=ray_origin(3)**2-cth2*origin2
4691 if (abs(aa)<1.d-14)
then
4694 disc=bb**2-4.d0*aa*cc
4695 if (disc>=
zero)
then
4696 root=sqrt(max(
zero,disc))
4702 do ix3=ixomin3,ixomax3+1
4704 if (abs(denom)>=1.d-14)
then
4705 numer=-phi_sin(ix3)*ray_origin(1)+phi_cos(ix3)*ray_origin(2)
4717 posmid=ray_origin+tmid*
vec_los
4720 if (dotp>=
zero)
then
4721 if (rc2<=rthick2) cycle
4723 if (r2<=rthick2) cycle
4728 mu=posmid(3)/sqrt(r2)
4733 phi=atan2(posmid(2),posmid(1))
4734 if (phi<phiface(ixomin3)-1.d-12) phi=phi+2.d0*
dpi
4735 if (phi>phiface(ixomax3+1)+1.d-12) phi=phi-2.d0*
dpi
4737 inside=ix1>0 .and. ix2>0 .and. ix3>0
4738 if (.not. inside) cycle
4739 ds_cm=(t1-t0)*unit_length
4740 if (si_unit) ds_cm=ds_cm*1.d2
4741 jds=max(
zero,source(ix^d))*ds_cm
4742 kds=max(
zero,opacity(ix^d))*ds_cm
4748 double precision,
intent(in) :: pos(1:3)
4749 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4750 phiface(ixOmin3:ixOmax3+1)
4751 integer,
intent(in) :: ixO^L
4752 integer,
intent(out) :: ix1,ix2,ix3
4753 logical,
intent(out) :: inside
4755 double precision :: r2,mu,phi
4765 phi=atan2(pos(2),pos(1))
4766 if (phi<phiface(ixomin3)-1.d-12) phi=phi+2.d0*dpi
4767 if (phi>phiface(ixomax3+1)+1.d-12) phi=phi-2.d0*dpi
4769 inside=ix1>0 .and. ix2>0 .and. ix3>0
4773 double precision,
intent(in) :: t,tNow,tExit,epsRay
4774 double precision,
intent(inout) :: tNext
4775 logical,
intent(inout) :: found
4777 if (t>tnow+epsray .and. t<=texit+epsray .and. t<tnext)
then
4784 double precision,
intent(in) :: t,theta_face_cos,ray_origin(1:3),tNow,tExit,epsRay
4785 double precision,
intent(inout) :: tNext
4786 logical,
intent(inout) :: found
4788 double precision :: pos(1:3),r2
4790 if (t<=tnow+epsray .or. t>texit+epsray .or. t>=tnext)
return
4793 if (r2<=zero)
return
4794 if (theta_face_cos>1.d-12 .and. pos(3)<-1.d-10)
return
4795 if (theta_face_cos<-1.d-12 .and. pos(3)>1.d-10)
return
4796 if (abs(pos(3)**2-theta_face_cos**2*r2)>1.d-6*max(one,r2))
return
4802 double precision,
intent(in) :: t,phi_face_sin,phi_face_cos,ray_origin(1:3),tNow,tExit,epsRay
4803 double precision,
intent(inout) :: tNext
4804 logical,
intent(inout) :: found
4806 double precision :: pos(1:3),radialDot
4808 if (t<=tnow+epsray .or. t>texit+epsray .or. t>=tnext)
return
4810 radialdot=phi_face_cos*pos(1)+phi_face_sin*pos(2)
4811 if (radialdot<-1.d-10)
return
4817 ix1,ix2,ix3,tNow,tExit,epsRay,tNext,found)
4818 double precision,
intent(in) :: ray_origin(1:3)
4819 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4820 phiface(ixOmin3:ixOmax3+1)
4821 double precision,
intent(in) :: phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
4822 integer,
intent(in) :: ixO^L,ix1,ix2,ix3
4823 double precision,
intent(in) :: tNow,tExit,epsRay
4824 double precision,
intent(out) :: tNext
4825 logical,
intent(out) :: found
4828 double precision :: dir2,origin2,odotd,aa,bb,cc,disc,root,cth2,denom,numer
4833 origin2=sum(ray_origin**2)
4839 cc=origin2-rface2(iface)
4840 disc=bb**2-4.d0*aa*cc
4841 if (disc>=zero)
then
4842 root=sqrt(max(zero,disc))
4849 cth2=theta_cos(iface)**2
4851 bb=2.d0*(ray_origin(3)*
vec_los(3)-cth2*odotd)
4852 cc=ray_origin(3)**2-cth2*origin2
4853 if (abs(aa)<1.d-14)
then
4855 ray_origin,tnow,texit,epsray,tnext,found)
4857 disc=bb**2-4.d0*aa*cc
4858 if (disc>=zero)
then
4859 root=sqrt(max(zero,disc))
4861 ray_origin,tnow,texit,epsray,tnext,found)
4863 ray_origin,tnow,texit,epsray,tnext,found)
4870 if (abs(denom)>=1.d-14)
then
4871 numer=-phi_sin(iface)*ray_origin(1)+phi_cos(iface)*ray_origin(2)
4873 tnow,texit,epsray,tnext,found)
4879 ray_origin,ximg1,ximg2,rface2,theta_cos,phiface,&
4881 t_enter,t_exit,segments,nseg,capacity,ok)
4884 integer,
intent(in) :: ixI^L,ixO^L,pixel_id
4885 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S)
4886 double precision,
intent(in) :: ray_origin(1:3),ximg1,ximg2
4887 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4888 phiface(ixOmin3:ixOmax3+1)
4889 double precision,
intent(in) :: phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
4890 double precision,
intent(in) :: t_enter,t_exit
4891 double precision,
allocatable,
intent(inout) :: segments(:,:)
4892 integer,
intent(inout) :: nseg,capacity
4893 logical,
intent(out) :: ok
4895 integer :: ix^D,nstep,maxSteps
4896 double precision :: tNow,tNext,tEnd,tMid,epsRay,ds_cm,jds,kds
4897 double precision :: pos(1:3),r2,dotp,rthick2,rc2
4898 logical :: inside,found
4901 if (t_exit<=t_enter)
return
4902 epsray=max(1.d-12,1.d-10*max(
one,abs(t_exit-t_enter)))
4903 pos=ray_origin+(t_enter+epsray)*
vec_los
4905 if (.not. inside)
then
4910 rthick2=(r_opt_thick*
const_rsun/unit_length)**2
4911 rc2=ximg1**2+ximg2**2
4914 maxsteps=8*((ixomax1-ixomin1+1)+(ixomax2-ixomin2+1)+(ixomax3-ixomin3+1)+3)
4916 do while (tnow<t_exit-epsray)
4918 ix1,ix2,ix3,tnow,t_exit,epsray,tnext,found)
4919 if (.not. found)
then
4923 tend=min(tnext,t_exit)
4925 tmid=
half*(tnow+tend)
4929 if (.not. ((dotp>=
zero .and. rc2<=rthick2) .or. &
4930 (dotp<
zero .and. r2<=rthick2)))
then
4931 ds_cm=(tend-tnow)*unit_length
4932 if (si_unit) ds_cm=ds_cm*1.d2
4933 jds=max(
zero,source(ix^d))*ds_cm
4934 kds=max(
zero,opacity(ix^d))*ds_cm
4939 if (tnow>=t_exit-epsray)
exit
4940 pos=ray_origin+(tnow+epsray)*
vec_los
4942 if (.not. inside)
then
4947 if (nstep>maxsteps)
then
4955 pixel_id,ray_origin,ximg1,ximg2,&
4956 rface,thetaface,phiface,rface2,&
4957 theta_cos,phi_sin,phi_cos,&
4958 segments,nseg,capacity,fallback)
4959 integer,
intent(in) :: ixI^L,ixO^L,pixel_id
4960 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S)
4961 double precision,
intent(in) :: ray_origin(1:3),ximg1,ximg2
4962 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4963 phiface(ixOmin3:ixOmax3+1)
4964 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4965 phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
4966 double precision,
allocatable,
intent(inout) :: segments(:,:)
4967 integer,
intent(inout) :: nseg,capacity
4968 logical,
intent(out) :: fallback
4970 integer :: nt,i,nsegStart,ix^D
4971 double precision :: tvals(12),posMid(1:3),t0,t1,tMid
4972 double precision :: dir2,origin2,odotd,aa,bb,cc,disc,root,cth2,denom,numer
4973 logical :: inside,ok
4979 origin2=sum(ray_origin**2)
4982 do ix1=ixomin1,ixomax1+1,ixomax1-ixomin1+1
4985 cc=origin2-rface2(ix1)
4986 disc=bb**2-4.d0*aa*cc
4987 if (disc>=zero)
then
4988 root=sqrt(max(zero,disc))
4993 do ix2=ixomin2,ixomax2+1,ixomax2-ixomin2+1
4994 cth2=theta_cos(ix2)**2
4996 bb=2.d0*(ray_origin(3)*
vec_los(3)-cth2*odotd)
4997 cc=ray_origin(3)**2-cth2*origin2
4998 if (abs(aa)<1.d-14)
then
5001 disc=bb**2-4.d0*aa*cc
5002 if (disc>=zero)
then
5003 root=sqrt(max(zero,disc))
5009 do ix3=ixomin3,ixomax3+1,ixomax3-ixomin3+1
5011 if (abs(denom)>=1.d-14)
then
5012 numer=-phi_sin(ix3)*ray_origin(1)+phi_cos(ix3)*ray_origin(2)
5025 posmid=ray_origin+tmid*
vec_los
5027 if (.not. inside) cycle
5029 ray_origin,ximg1,ximg2,rface2,theta_cos,phiface,phi_sin,phi_cos,&
5030 t0,t1,segments,nseg,capacity,ok)
5035 pixel_id,ray_origin,ximg1,ximg2,rface,thetaface,phiface,rface2,&
5036 theta_cos,phi_sin,phi_cos,segments,nseg,capacity)
5046 integer,
intent(in) :: numXI1,numXI2
5047 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2),dxI
5049 double precision,
intent(inout) :: EM(numXI1,numXI2)
5051 integer :: ixO^L,ixI^L,ix^D
5052 integer :: iigrid,igrid,ixP1,ixP2,ixPmin1,ixPmax1,ixPmin2,ixPmax2
5053 integer :: iseg,nseg,capacity,sphDdaFallbackLocal,sphDdaFallbackGlobal
5054 double precision :: ray_origin(1:3),profile_local(3),profile_global(3)
5055 double precision,
allocatable :: source(:^D&)
5056 double precision,
allocatable :: rface(:),thetaface(:),phiface(:)
5057 double precision,
allocatable :: rface2(:),theta_cos(:),phi_sin(:),phi_cos(:)
5058 double precision,
allocatable :: segments(:,:)
5059 logical :: has_pixels,ddaFallback
5062 sphddafallbacklocal=0
5063 do iigrid=1,igridstail; igrid=igrids(iigrid);
5064 ^d&ixomin^d=ixmlo^d\
5065 ^d&ixomax^d=ixmhi^d\
5066 ^d&iximin^d=
ixglo^d\
5067 ^d&iximax^d=
ixghi^d\
5069 allocate(source(ixi^s))
5074 allocate(rface2(ixomin1:ixomax1+1),theta_cos(ixomin2:ixomax2+1),&
5075 phi_sin(ixomin3:ixomax3+1),phi_cos(ixomin3:ixomax3+1))
5077 theta_cos=cos(thetaface)
5078 phi_sin=sin(phiface)
5079 phi_cos=cos(phiface)
5082 ixpmin1,ixpmax1,ixpmin2,ixpmax2,has_pixels)
5083 if (has_pixels)
then
5085 do ixp1=ixpmin1,ixpmax1
5086 do ixp2=ixpmin2,ixpmax2
5088 profile_local(1)=profile_local(1)+one
5092 1,ray_origin,xi1(ixp1),xi2(ixp2),rface,thetaface,phiface,&
5093 rface2,theta_cos,phi_sin,phi_cos,segments,nseg,capacity,ddafallback)
5094 if (ddafallback) sphddafallbacklocal=sphddafallbacklocal+1
5096 em(ixp1,ixp2)=em(ixp1,ixp2)+segments(3,iseg)
5100 rface,thetaface,phiface,em(ixp1,ixp2))
5104 profile_local(2)=profile_local(2)+dble((ixpmax1-ixpmin1+1)*(ixpmax2-ixpmin2+1))
5106 if (
allocated(segments))
deallocate(segments)
5107 profile_local(3)=profile_local(3)+one
5108 if (
allocated(rface2))
deallocate(rface2)
5109 if (
allocated(theta_cos))
deallocate(theta_cos)
5110 if (
allocated(phi_sin))
deallocate(phi_sin)
5111 if (
allocated(phi_cos))
deallocate(phi_cos)
5112 deallocate(source,rface,thetaface,phiface)
5114 call mpi_allreduce(profile_local,profile_global,3,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5115 call mpi_allreduce(sphddafallbacklocal,sphddafallbackglobal,1,mpi_integer,mpi_sum,
icomm,
ierrmpi)
5118 write(*,
'(a,3(es12.5,1x))')
' sph_dda thin profile rays pixels blocks: ',profile_global
5119 write(*,
'(a,i0)')
' sph_dda thin fallback rays: ',sphddafallbackglobal
5121 write(*,
'(a,3(es12.5,1x))')
' sph_intersection thin profile rays pixels blocks: ',profile_global
5129 integer,
intent(in) :: numXI1,numXI2
5130 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2),dxI
5132 double precision,
intent(out) :: EUV(numXI1,numXI2),Tau(numXI1,numXI2),EUVthin(numXI1,numXI2)
5134 integer,
parameter :: nSegVars=5
5135 integer :: ixO^L,ixI^L,ix^D
5136 integer :: iigrid,igrid,ixP1,ixP2,ipix,ipixStart,ipixEnd,nPixBatch,pixel_id
5137 integer :: nseg,capacity,totalCount,totalSeg,ipe,is,iseg,nidx,owner,isegDest,nsegBefore
5138 integer :: ixGlobal,iyGlobal,ixPmin1,ixPmax1,ixPmin2,ixPmax2,iFirst,iLast,iLocal
5139 integer :: nPixBatchTarget
5140 integer :: maxSegBatchTarget,maxSegCommTarget,maxNsegBatch,nPixTotal
5141 integer :: maxOwnerSegCount,maxOwnerSegCountLocal,segOffset,recvFill,totalRoundCount,totalRoundSeg
5142 integer :: sphDdaFallbackLocal,sphDdaFallbackGlobal
5143 integer,
allocatable :: sendCounts(:),recvCounts(:),sendDispls(:),recvDispls(:)
5144 integer,
allocatable :: roundSendCounts(:),roundRecvCounts(:)
5145 integer,
allocatable :: roundSendDispls(:),roundRecvDispls(:)
5146 integer,
allocatable :: ownerSegCounts(:),ownerOffsets(:),idx(:)
5147 integer,
allocatable :: bucketCounts(:),bucketOffsets(:),bucketFill(:)
5148 double precision :: ray_origin(1:3),atten
5149 double precision :: profile_local(5),profile_global(5),profile_batch(5)
5150 double precision :: phys_max_local(2),phys_max_global(2)
5151 double precision :: phys_sum_local(2),phys_sum_global(2),phys_sum_batch(2)
5152 double precision,
allocatable :: segments(:,:),segments_send(:,:),segments_recv(:,:)
5153 double precision,
allocatable :: segments_recv_round(:,:)
5154 double precision,
allocatable :: image_reduce(:,:)
5155 logical :: has_pixels,batchAccepted,batchReduced,ddaFallback
5164 sphddafallbacklocal=0
5165 allocate(sendcounts(0:
npe-1),recvcounts(0:
npe-1),senddispls(0:
npe-1),recvdispls(0:
npe-1))
5166 allocate(roundsendcounts(0:
npe-1),roundrecvcounts(0:
npe-1))
5167 allocate(roundsenddispls(0:
npe-1),roundrecvdispls(0:
npe-1))
5168 allocate(ownersegcounts(0:
npe-1),owneroffsets(0:
npe-1))
5169 allocate(cache(igridstail))
5171 allocate(bucketcounts(npixbatchtarget),bucketoffsets(npixbatchtarget+1),&
5172 bucketfill(npixbatchtarget))
5174 do iigrid=1,igridstail; igrid=igrids(iigrid);
5175 ^d&ixomin^d=ixmlo^d\
5176 ^d&ixomax^d=ixmhi^d\
5177 ^d&iximin^d=
ixglo^d\
5178 ^d&iximax^d=
ixghi^d\
5180 cache(iigrid)%igrid=igrid
5181 allocate(cache(iigrid)%source(ixi^s),cache(iigrid)%opacity(ixi^s))
5182 cache(iigrid)%source=zero
5183 cache(iigrid)%opacity=zero
5184 call get_euv(
wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,cache(iigrid)%source)
5187 phys_max_local(1)=max(phys_max_local(1),maxval(cache(iigrid)%source(ixo^s)))
5188 phys_max_local(2)=max(phys_max_local(2),maxval(cache(iigrid)%opacity(ixo^s)))
5190 cache(iigrid)%rface,cache(iigrid)%thetaface,&
5191 cache(iigrid)%phiface)
5192 allocate(cache(iigrid)%rface2(ixomin1:ixomax1+1),&
5193 cache(iigrid)%theta_cos(ixomin2:ixomax2+1),&
5194 cache(iigrid)%phi_sin(ixomin3:ixomax3+1),&
5195 cache(iigrid)%phi_cos(ixomin3:ixomax3+1))
5196 cache(iigrid)%rface2=cache(iigrid)%rface**2
5197 cache(iigrid)%theta_cos=cos(cache(iigrid)%thetaface)
5198 cache(iigrid)%phi_sin=sin(cache(iigrid)%phiface)
5199 cache(iigrid)%phi_cos=cos(cache(iigrid)%phiface)
5201 cache(iigrid)%phiface,ixo^l,numxi1,numxi2,xi1,xi2,dxi,&
5202 cache(iigrid)%ixPmin1,cache(iigrid)%ixPmax1,&
5203 cache(iigrid)%ixPmin2,cache(iigrid)%ixPmax2,cache(iigrid)%has_pixels)
5206 npixtotal=numxi1*numxi2
5208 do while (ipixstart<=npixtotal)
5209 ipixend=min(numxi1*numxi2,ipixstart+npixbatchtarget-1)
5210 npixbatch=ipixend-ipixstart+1
5211 batchaccepted=.false.
5212 batchreduced=.false.
5214 do while (.not. batchaccepted)
5220 do iigrid=1,igridstail; igrid=igrids(iigrid);
5221 ^d&ixomin^d=ixmlo^d\
5222 ^d&ixomax^d=ixmhi^d\
5223 ^d&iximin^d=
ixglo^d\
5224 ^d&iximax^d=
ixghi^d\
5226 ixpmin1=cache(iigrid)%ixPmin1
5227 ixpmax1=cache(iigrid)%ixPmax1
5228 ixpmin2=cache(iigrid)%ixPmin2
5229 ixpmax2=cache(iigrid)%ixPmax2
5230 has_pixels=cache(iigrid)%has_pixels
5231 if (.not. has_pixels) cycle
5233 do ixp2=ixpmin2,ixpmax2
5234 ifirst=max(ipixstart,(ixp2-1)*numxi1+ixpmin1)
5235 ilast=min(ipixend,(ixp2-1)*numxi1+ixpmax1)
5236 if (ifirst>ilast) cycle
5237 do ipix=ifirst,ilast
5238 ixp1=1+mod(ipix-1,numxi1)
5241 profile_batch(1)=profile_batch(1)+one
5245 cache(iigrid)%opacity,pixel_id,ray_origin,xi1(ixp1),xi2(ixp2),&
5246 cache(iigrid)%rface,cache(iigrid)%thetaface,cache(iigrid)%phiface,&
5247 cache(iigrid)%rface2,cache(iigrid)%theta_cos,&
5248 cache(iigrid)%phi_sin,cache(iigrid)%phi_cos,&
5249 segments,nseg,capacity,ddafallback)
5250 if (ddafallback) sphddafallbacklocal=sphddafallbacklocal+1
5253 cache(iigrid)%opacity,pixel_id,ray_origin,xi1(ixp1),xi2(ixp2),&
5254 cache(iigrid)%rface,cache(iigrid)%thetaface,cache(iigrid)%phiface,&
5255 cache(iigrid)%rface2,cache(iigrid)%theta_cos,&
5256 cache(iigrid)%phi_sin,cache(iigrid)%phi_cos,&
5257 segments,nseg,capacity)
5259 if (nseg>nsegbefore) profile_batch(2)=profile_batch(2)+one
5260 profile_batch(3)=profile_batch(3)+dble(nseg-nsegbefore)
5261 do iseg=nsegbefore+1,nseg
5262 phys_sum_batch(1)=phys_sum_batch(1)+segments(3,iseg)
5263 phys_sum_batch(2)=phys_sum_batch(2)+segments(4,iseg)
5269 call mpi_allreduce(nseg,maxnsegbatch,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
5270 if (maxnsegbatch>maxsegbatchtarget .and. npixbatch>1)
then
5271 npixbatch=max(1,npixbatch/2)
5272 ipixend=ipixstart+npixbatch-1
5273 if (
allocated(segments))
deallocate(segments)
5276 batchaccepted=.true.
5280 profile_local=profile_local+profile_batch
5281 phys_sum_local=phys_sum_local+phys_sum_batch
5284 write(*,
'(a,3(i0,1x))')
' sph_dda thick adaptive batch: ',&
5285 ipixstart,ipixend,maxnsegbatch
5287 write(*,
'(a,3(i0,1x))')
' sph_intersection thick adaptive batch: ',&
5288 ipixstart,ipixend,maxnsegbatch
5292 if (.not.
allocated(segments))
then
5294 allocate(segments(nsegvars,capacity))
5299 ownersegcounts(owner)=ownersegcounts(owner)+1
5301 sendcounts=nsegvars*ownersegcounts
5304 senddispls(ipe)=senddispls(ipe-1)+sendcounts(ipe-1)
5307 allocate(segments_send(nsegvars,max(1,nseg)))
5311 isegdest=senddispls(owner)/nsegvars+owneroffsets(owner)+1
5312 segments_send(:,isegdest)=segments(:,is)
5313 owneroffsets(owner)=owneroffsets(owner)+1
5316 call mpi_alltoall(sendcounts,1,mpi_integer,recvcounts,1,mpi_integer,
icomm,
ierrmpi)
5319 recvdispls(ipe)=recvdispls(ipe-1)+recvcounts(ipe-1)
5321 totalcount=sum(recvcounts)
5322 totalseg=totalcount/nsegvars
5323 profile_local(4)=profile_local(4)+dble(totalcount)
5324 allocate(segments_recv(nsegvars,max(1,totalseg)))
5327 maxownersegcountlocal=maxval(ownersegcounts)
5328 call mpi_allreduce(maxownersegcountlocal,maxownersegcount,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
5329 do segoffset=0,maxownersegcount-1,maxsegcommtarget
5331 roundsenddispls=senddispls
5333 if (ownersegcounts(ipe)>segoffset)
then
5334 roundsendcounts(ipe)=nsegvars*min(maxsegcommtarget,ownersegcounts(ipe)-segoffset)
5335 roundsenddispls(ipe)=senddispls(ipe)+nsegvars*segoffset
5339 call mpi_alltoall(roundsendcounts,1,mpi_integer,roundrecvcounts,1,mpi_integer,
icomm,
ierrmpi)
5340 roundrecvdispls(0)=0
5342 roundrecvdispls(ipe)=roundrecvdispls(ipe-1)+roundrecvcounts(ipe-1)
5344 totalroundcount=sum(roundrecvcounts)
5345 totalroundseg=totalroundcount/nsegvars
5346 allocate(segments_recv_round(nsegvars,max(1,totalroundseg)))
5348 call mpi_alltoallv(segments_send,roundsendcounts,roundsenddispls,mpi_double_precision,&
5349 segments_recv_round,roundrecvcounts,roundrecvdispls,&
5352 if (totalroundseg>0)
then
5353 segments_recv(:,recvfill+1:recvfill+totalroundseg)=segments_recv_round(:,1:totalroundseg)
5354 recvfill=recvfill+totalroundseg
5356 deallocate(segments_recv_round)
5359 if (recvfill/=totalseg)
call mpistop(
"ray-segment receive mismatch")
5361 if (totalseg>0)
then
5362 allocate(idx(totalseg))
5363 bucketcounts(1:npixbatch)=0
5366 ipix=nint(segments_recv(1,is))
5368 ilocal=ipix-ipixstart+1
5369 bucketcounts(ilocal)=bucketcounts(ilocal)+1
5375 do ilocal=1,npixbatch
5376 bucketoffsets(ilocal+1)=bucketoffsets(ilocal)+bucketcounts(ilocal)
5378 bucketfill(1:npixbatch)=bucketoffsets(1:npixbatch)
5381 ipix=nint(segments_recv(1,is))
5383 ilocal=ipix-ipixstart+1
5384 idx(bucketfill(ilocal))=is
5385 bucketfill(ilocal)=bucketfill(ilocal)+1
5390 do ipix=ipixstart,ipixend
5392 ilocal=ipix-ipixstart+1
5393 nidx=bucketcounts(ilocal)
5395 profile_local(5)=profile_local(5)+dble(nidx)*dble(nidx)
5397 idx(bucketoffsets(ilocal):bucketoffsets(ilocal+1)-1),nidx)
5398 ixglobal=1+mod(ipix-1,numxi1)
5399 iyglobal=1+(ipix-1)/numxi1
5400 do iseg=bucketoffsets(ilocal),bucketoffsets(ilocal+1)-1
5402 euvthin(ixglobal,iyglobal)=euvthin(ixglobal,iyglobal)+segments_recv(3,is)
5404 euv(ixglobal,iyglobal)=euv(ixglobal,iyglobal)+atten*segments_recv(3,is)
5405 tau(ixglobal,iyglobal)=tau(ixglobal,iyglobal)+max(zero,segments_recv(4,is))
5412 deallocate(segments_send,segments_recv)
5413 if (
allocated(segments))
deallocate(segments)
5417 do iigrid=1,igridstail
5418 if (
allocated(cache(iigrid)%source))
deallocate(cache(iigrid)%source)
5419 if (
allocated(cache(iigrid)%opacity))
deallocate(cache(iigrid)%opacity)
5420 if (
allocated(cache(iigrid)%rface))
deallocate(cache(iigrid)%rface)
5421 if (
allocated(cache(iigrid)%thetaface))
deallocate(cache(iigrid)%thetaface)
5422 if (
allocated(cache(iigrid)%phiface))
deallocate(cache(iigrid)%phiface)
5423 if (
allocated(cache(iigrid)%rface2))
deallocate(cache(iigrid)%rface2)
5424 if (
allocated(cache(iigrid)%theta_cos))
deallocate(cache(iigrid)%theta_cos)
5425 if (
allocated(cache(iigrid)%phi_sin))
deallocate(cache(iigrid)%phi_sin)
5426 if (
allocated(cache(iigrid)%phi_cos))
deallocate(cache(iigrid)%phi_cos)
5429 deallocate(sendcounts,recvcounts,senddispls,recvdispls,roundsendcounts,roundrecvcounts,&
5430 roundsenddispls,roundrecvdispls,ownersegcounts,owneroffsets,bucketcounts,&
5431 bucketoffsets,bucketfill)
5432 allocate(image_reduce(numxi1,numxi2))
5433 call mpi_allreduce(euv,image_reduce,numxi1*numxi2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5435 call mpi_allreduce(tau,image_reduce,numxi1*numxi2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5437 call mpi_allreduce(euvthin,image_reduce,numxi1*numxi2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5438 euvthin=image_reduce
5439 deallocate(image_reduce)
5440 call mpi_allreduce(profile_local,profile_global,5,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5441 call mpi_allreduce(phys_max_local,phys_max_global,2,mpi_double_precision,mpi_max,
icomm,
ierrmpi)
5442 call mpi_allreduce(phys_sum_local,phys_sum_global,2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5443 call mpi_allreduce(sphddafallbacklocal,sphddafallbackglobal,1,mpi_integer,mpi_sum,
icomm,
ierrmpi)
5446 write(*,
'(a,5(es12.5,1x))')
' sph_dda thick profile: ',profile_global
5447 write(*,
'(a,i0)')
' sph_dda thick fallback rays: ',sphddafallbackglobal
5448 write(*,
'(a,4(es12.5,1x))')
' sph_dda thick physics maxj maxk sumjds sumkds: ',&
5449 phys_max_global(1),phys_max_global(2),phys_sum_global(1),phys_sum_global(2)
5451 write(*,
'(a,5(es12.5,1x))')
' sph_intersection thick profile: ',profile_global
5452 write(*,
'(a,4(es12.5,1x))')
' sph_intersection thick physics maxj maxk sumjds sumkds: ',&
5453 phys_max_global(1),phys_max_global(2),phys_sum_global(1),phys_sum_global(2)
5463 double precision,
intent(out) :: xImin1,xImax1,xImin2,xImax2
5465 integer,
parameter :: nsample=5
5466 integer :: iigrid,igrid,ir,it,ip
5467 integer :: ixI^L,ixO^L
5468 double precision,
allocatable :: rface(:),thetaface(:),phiface(:)
5469 double precision :: local_min1,local_max1,local_min2,local_max2
5470 double precision :: sph(1:3),xcent(1:2),wr,wt,wp
5472 local_min1=huge(one)
5473 local_max1=-huge(one)
5474 local_min2=huge(one)
5475 local_max2=-huge(one)
5477 do iigrid=1,igridstail
5478 igrid=igrids(iigrid)
5479 ^d&ixomin^d=ixmlo^d\
5480 ^d&ixomax^d=ixmhi^d\
5481 ^d&iximin^d=ixglo^d\
5482 ^d&iximax^d=ixghi^d\
5485 rface,thetaface,phiface)
5487 wr=dble(ir)/dble(nsample-1)
5488 sph(1)=(one-wr)*rface(ixomin1)+wr*rface(ixomax1+1)
5490 wt=dble(it)/dble(nsample-1)
5491 sph(2)=(one-wt)*thetaface(ixomin2)+wt*thetaface(ixomax2+1)
5493 wp=dble(ip)/dble(nsample-1)
5494 if (ir/=0 .and. ir/=nsample-1 .and. it/=0 .and. it/=nsample-1 .and. &
5495 ip/=0 .and. ip/=nsample-1) cycle
5496 sph(3)=(one-wp)*phiface(ixomin3)+wp*phiface(ixomax3+1)
5498 local_min1=min(local_min1,xcent(1))
5499 local_max1=max(local_max1,xcent(1))
5500 local_min2=min(local_min2,xcent(2))
5501 local_max2=max(local_max2,xcent(2))
5505 deallocate(rface,thetaface,phiface)
5508 call mpi_allreduce(local_min1,ximin1,1,mpi_double_precision,mpi_min,icomm,ierrmpi)
5509 call mpi_allreduce(local_max1,ximax1,1,mpi_double_precision,mpi_max,icomm,ierrmpi)
5510 call mpi_allreduce(local_min2,ximin2,1,mpi_double_precision,mpi_min,icomm,ierrmpi)
5511 call mpi_allreduce(local_max2,ximax2,1,mpi_double_precision,mpi_max,icomm,ierrmpi)
5512 if (ximin1>0.5d0*huge(one) .or. ximax1<-0.5d0*huge(one) .or. &
5513 ximin2>0.5d0*huge(one) .or. ximax2<-0.5d0*huge(one))
then
5514 call mpistop(
"sph_intersection could not determine image bounds")
5519 double precision,
intent(out) :: dxI
5521 integer :: iigrid,igrid,ixI^L,ixO^L,ix^D
5522 double precision :: local_min,global_min,dr,ds_theta,ds_phi,rval,theta
5525 do iigrid=1,igridstail
5526 igrid=igrids(iigrid)
5527 ^d&ixomin^d=ixmlo^d\
5528 ^d&ixomax^d=ixmhi^d\
5529 ^d&iximin^d=ixglo^d\
5530 ^d&iximax^d=ixghi^d\
5532 do ix1=ixomin1,ixomax1
5533 do ix2=ixomin2,ixomax2
5534 do ix3=ixomin3,ixomax3
5535 rval=max(smalldouble,ps(igrid)%x(ix^d,1))
5536 theta=ps(igrid)%x(ix^d,2)
5537 dr=ps(igrid)%dx(ix^d,1)
5538 ds_theta=rval*ps(igrid)%dx(ix^d,2)
5539 ds_phi=rval*max(smalldouble,sin(theta))*ps(igrid)%dx(ix^d,3)
5540 local_min=min(local_min,dr,ds_theta,ds_phi)
5546 call mpi_allreduce(local_min,global_min,1,mpi_double_precision,mpi_min,icomm,ierrmpi)
5547 if (global_min<=zero .or. global_min>half*huge(one))
then
5548 call mpistop(
"sph_intersection could not determine dat-resolution image spacing")
5559 integer,
intent(in) :: qunit
5561 character(20),
intent(in) :: datatype
5563 integer :: ix^D,numXI1,numXI2,numWI
5564 double precision :: xImin1,xImax1,xImin2,xImax2,xIcent1,xIcent2,dxI
5565 double precision,
allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:)
5566 double precision,
allocatable :: wI(:,:,:),wIs(:,:,:),EM(:,:),Dpl(:,:),Tau(:,:),EMthin(:,:),WLB(:,:,:)
5567 double precision :: vec_temp1(1:3),vec_temp2(1:3)
5568 double precision :: vec_z(1:3),vec_cor(1:3),xI_cor(1:2)
5569 double precision :: res,LOS_psi,r_max,r_loc
5572 character (30) :: ion
5573 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
5574 double precision :: arcsec,RHESSI_rsl,LASCO_rsl,pixel,R_occult,smallflux
5575 integer :: iigrid,igrid,i,j,numSI,iw
5576 logical :: emit,ray_image_global,has_thick_output
5578 if (coordinate==spherical)
then
5586 if (coordinate==spherical)
then
5590 ximin1=-abs(xprobmax1)
5591 ximin2=-abs(xprobmax1)
5592 ximax1=abs(xprobmax1)
5593 ximax2=abs(xprobmax1)
5598 if (ix1==1) vec_cor(1)=xprobmin1
5599 if (ix1==2) vec_cor(1)=xprobmax1
5601 if (ix2==1) vec_cor(2)=xprobmin2
5602 if (ix2==2) vec_cor(2)=xprobmax2
5604 if (ix3==1) vec_cor(3)=xprobmin3
5605 if (ix3==2) vec_cor(3)=xprobmax3
5608 r_loc=r_loc+(vec_cor(2)-
x_origin(2))**2
5609 r_loc=r_loc+(vec_cor(3)-
x_origin(3))**2
5611 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
5614 r_max=max(r_max,r_loc)
5618 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
5624 ximin1=min(ximin1,xi_cor(1))
5625 ximax1=max(ximax1,xi_cor(1))
5626 ximin2=min(ximin2,xi_cor(2))
5627 ximax2=max(ximax2,xi_cor(2))
5640 xicent1=(ximin1+ximax1)/2.d0
5641 xicent2=(ximin2+ximax2)/2.d0
5649 if (datatype==
'image_euv')
then
5653 else if (datatype==
'image_sxr')
then
5655 dxi=rhessi_rsl*arcsec
5657 else if (datatype==
'image_whitelight')
then
5668 call mpistop(
'Whitelight synthesis: instrument is not supported!')
5670 dxi=lasco_rsl*arcsec
5675 numxi1=8*ceiling((ximax1-xicent1)/dxi/8.d0)
5676 ximin1=xicent1-numxi1*dxi
5677 ximax1=xicent1+numxi1*dxi
5679 numxi2=8*ceiling((ximax2-xicent2)/dxi/8.d0)
5680 ximin2=xicent2-numxi2*dxi
5681 ximax2=xicent2+numxi2*dxi
5683 allocate(xi1(numxi1),xi2(numxi2),dxi1(numxi1),dxi2(numxi2))
5685 xi1(ix1)=ximin1+dxi*(ix1-
half)
5689 xi2(ix2)=ximin2+dxi*(ix2-
half)
5694 if (datatype==
'image_euv' .or. datatype==
'image_sxr')
then
5695 has_thick_output=datatype==
'image_euv' .and. trim(
radiation_transfer)==
'thick' .and. &
5698 if (datatype==
'image_euv')
then
5703 allocate(wi(numxi1,numxi2,numwi),wis(numxi1,numxi2,numwi),em(numxi1,numxi2))
5707 ray_image_global=.false.
5708 if (has_thick_output)
then
5709 allocate(tau(numxi1,numxi2),emthin(numxi1,numxi2))
5713 if (coordinate==cartesian .and. datatype==
'image_euv' .and. &
5715 ray_image_global=.true.
5716 allocate(dpl(numxi1,numxi2))
5724 else if (coordinate==cartesian)
then
5725 do iigrid=1,igridstail; igrid=igrids(iigrid);
5728 else if (trim(
ray_method_active) ==
'spherical' .and. datatype ==
'image_euv')
then
5730 ray_image_global=.true.
5736 do iigrid=1,igridstail; igrid=igrids(iigrid);
5740 if (ray_image_global)
then
5741 if (has_thick_output)
then
5743 has_thick_output,tau=tau,euvthin=emthin,&
5744 cap_absorption=.true.)
5751 if (em(ix1,ix2)>smallflux) wis(ix1,ix2,1)=em(ix1,ix2)
5755 if (.not. ray_image_global)
then
5756 numsi=numxi1*numxi2*numwi
5757 call mpi_allreduce(wis,wi,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5765 call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
5766 if (
allocated(tau))
deallocate(tau)
5767 if (
allocated(emthin))
deallocate(emthin)
5768 deallocate(wi,wis,em)
5769 else if (datatype==
'image_whitelight')
then
5771 allocate(wi(numxi1,numxi2,numwi),wis(numxi1,numxi2,numwi),wlb(numxi1,numxi2,numwi))
5775 if (coordinate==spherical)
then
5776 do iigrid=1,igridstail; igrid=igrids(iigrid);
5782 if (wlb(ix1,ix2,1)>smallflux)
then
5783 wis(ix1,ix2,1)=wlb(ix1,ix2,1)
5784 wis(ix1,ix2,2)=wlb(ix1,ix2,2)
5788 numsi=numxi1*numxi2*numwi
5789 call mpi_allreduce(wis,wi,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5796 call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
5797 deallocate(wi,wis,wlb)
5800 deallocate(xi1,xi2,dxi1,dxi2)
5805 integer,
intent(in) :: igrid,numXI1,numXI2
5806 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
5807 double precision,
intent(in) :: dxI
5809 character(20),
intent(in) :: datatype
5810 double precision,
intent(inout) :: EM(numXI1,numXI2)
5812 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
5813 double precision :: xb^L,xd^D
5814 double precision,
allocatable :: flux(:^D&),opacity(:^D&)
5815 double precision :: res
5816 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
5817 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC
5818 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2)
5821 double precision :: logTe
5822 character (30) :: ion
5823 double precision :: lineCent
5824 double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
5825 double precision :: arcsec,pixel,RHESSI_rsl,area_1AU
5826 double precision :: aa,bb
5828 ^d&ixomin^d=ixmlo^d\
5829 ^d&ixomax^d=ixmhi^d\
5830 ^d&iximin^d=ixglo^d\
5831 ^d&iximax^d=ixghi^d\
5832 ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
5833 ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
5836 arcsec=7.25d5/unit_length
5838 arcsec=7.25d7/unit_length
5841 allocate(flux(ixi^s),opacity(ixi^s))
5842 if (datatype==
'image_euv')
then
5843 if (trim(emission_model)==
'pseudo_current')
then
5845 else if (trim(emission_model)==
'radio_ff')
then
5849 call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
5850 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
5852 call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
5853 pixel=spacersl*arcsec
5854 sigma0=sigma_psf*pixel
5855 else if (datatype==
'image_sxr')
then
5857 call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr)
5858 rhessi_rsl=2.3d0/instrument_resolution_factor
5860 pixel=rhessi_rsl*arcsec
5861 sigma0=sigma_psf*pixel
5866 {
do ix^d=ixomin^d,ixomax^d\}
5868 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi1(^d))/(dxi/2.d0)));
5869 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi2(^d))/(dxi/2.d0)));
5870 ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
5871 if (datatype==
'image_euv')
then
5873 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length*1.d2/dxi/dxi
5875 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length/dxi/dxi
5877 else if (datatype==
'image_sxr')
then
5879 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length**3/area_1au
5881 if (fluxsubc>smalldouble)
then
5883 {
do isubc^d=1,nsubc^d\}
5884 ^d&xsubc(^d)=ps(igrid)%x(ix^dd,^d)-half*ps(igrid)%dx(ix^dd,^d)+(isubc^d-half)*dxsubc^d;
5888 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
5889 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
5890 ixpmin1=max(1,ixp1-3)
5891 ixpmax1=min(ixp1+3,numxi1)
5892 ixpmin2=max(1,ixp2-3)
5893 ixpmax2=min(ixp2+3,numxi2)
5894 do ixp1=ixpmin1,ixpmax1
5895 do ixp2=ixpmin2,ixpmax2
5896 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
5897 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
5898 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
5899 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
5900 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
5901 em(ixp1,ixp2)=em(ixp1,ixp2)+fluxsubc*factor
5908 deallocate(flux,opacity)
5912 integer,
intent(in) :: igrid,numXI1,numXI2
5913 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
5914 double precision,
intent(in) :: dxI
5916 character(20),
intent(in) :: datatype
5917 double precision,
intent(inout) :: EM(numXI1,numXI2)
5919 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
5920 double precision,
allocatable :: flux(:^D&),Ne(:^D&),opacity(:^D&)
5921 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
5922 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC,RsubC
5923 double precision :: TBsubC,PBsubC
5924 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2),xSubC_car(1:3)
5925 double precision :: R_thick,dotp,dvolume,R_occult,Rc
5926 double precision :: dxl(1:3),x_sph(1:3),dx_sph(1:3)
5927 double precision :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
5928 logical :: sun_back_side,emit
5931 double precision :: logTe
5932 character (30) :: ion
5933 double precision :: lineCent
5934 double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
5935 double precision :: RHESSI_rsl,area_1AU,arcsec,pixel
5937 ^d&ixomin^d=ixmlo^d;
5938 ^d&ixomax^d=ixmhi^d;
5939 ^d&iximin^d=ixglo^d;
5940 ^d&iximax^d=ixghi^d;
5943 arcsec=7.25d5/unit_length
5945 arcsec=7.25d7/unit_length
5948 allocate(flux(ixi^s),opacity(ixi^s))
5949 if (datatype==
'image_euv')
then
5950 if (trim(emission_model)==
'pseudo_current')
then
5952 else if (trim(emission_model)==
'radio_ff')
then
5956 call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
5957 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
5959 call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
5960 pixel=spacersl*arcsec
5961 sigma0=sigma_psf*pixel
5962 else if (datatype==
'image_sxr')
then
5964 call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr)
5965 rhessi_rsl=2.3d0/instrument_resolution_factor
5967 pixel=rhessi_rsl*arcsec
5968 sigma0=sigma_psf*pixel
5973 r_thick=r_opt_thick*const_rsun/unit_length
5974 {
do ix^d=ixomin^d,ixomax^d\}
5975 x_sph(1:3)=ps(igrid)%x(ix^d,1:3)
5976 dx_sph(1:3)=ps(igrid)%dx(ix^d,1:3)
5978 dxl(2)=x_sph(1)*dx_sph(2)
5979 dxl(3)=x_sph(1)*dsin(x_sph(2))*dx_sph(3)
5984 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
5986 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
5988 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
5990 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
5992 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
5994 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
5999 xsubc(1)=x_sph(1)-half*dx_sph(1)+(isubc1-half)*dx_sph(1)/nsubc1
6001 dxsubc1=dx_sph(1)/nsubc1
6004 xsubc(2)=x_sph(2)-half*dx_sph(2)+(isubc2-half)*dx_sph(2)/nsubc2
6005 dxsubc2=xsubc(1)*dx_sph(2)/nsubc2
6006 dxsubc3=xsubc(1)*dsin(xsubc(2))*dx_sph(3)/nsubc3
6007 dvolume=dxsubc1*dxsubc2*dxsubc3
6008 if (datatype==
'image_euv')
then
6010 fluxsubc=flux(ix^d)*dvolume*unit_length*1.d2/dxi/dxi
6012 fluxsubc=flux(ix^d)*dvolume*unit_length/dxi/dxi
6014 else if (datatype==
'image_sxr')
then
6016 fluxsubc=flux(ix^d)*dvolume*unit_length**3/area_1au
6019 if (fluxsubc>smalldouble)
then
6022 xsubc(3)=x_sph(3)-half*dx_sph(3)+(isubc3-half)*dx_sph(3)/nsubc3
6024 rc=dsqrt(xcent(1)**2+xcent(2)**2)
6029 sun_back_side=.true.
6030 if (dotp<0.d0) sun_back_side=.false.
6032 if (sun_back_side)
then
6034 if (rc>r_thick) emit=.true.
6037 if (xsubc(1)<=r_thick) emit=.false.
6043 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
6044 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
6045 ixpmin1=max(1,ixp1-3)
6046 ixpmax1=min(ixp1+3,numxi1)
6047 ixpmin2=max(1,ixp2-3)
6048 ixpmax2=min(ixp2+3,numxi2)
6049 do ixp1=ixpmin1,ixpmax1
6050 do ixp2=ixpmin2,ixpmax2
6051 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6052 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6053 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6054 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6055 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
6056 em(ixp1,ixp2)=em(ixp1,ixp2)+fluxsubc*factor
6066 deallocate(flux,opacity)
6073 integer,
intent(in) :: igrid,numXI1,numXI2,numWI
6074 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
6075 double precision,
intent(in) :: dxI
6077 character(20),
intent(in) :: datatype
6078 double precision,
intent(inout) :: WLB(numXI1,numXI2,numWI)
6080 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
6081 double precision,
allocatable :: flux(:^D&),Ne(:^D&)
6082 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
6083 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC,RsubC
6084 double precision :: sigma_PSF,sigma0,arcsec,pixel,LASCO_rsl
6085 double precision :: A,B,C,D,Rc,Ne0,TBsubC,PBsubC,factor
6086 double precision :: R_thick,dotp,dvolume,R_occult
6087 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2),xSubC_car(1:3)
6088 double precision :: dxl(1:3),x_sph(1:3),dx_sph(1:3)
6089 double precision :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
6092 ^d&ixomin^d=ixmlo^d;
6093 ^d&ixomax^d=ixmhi^d;
6094 ^d&iximin^d=ixglo^d;
6095 ^d&iximax^d=ixghi^d;
6098 arcsec=7.25d5/unit_length
6100 arcsec=7.25d7/unit_length
6104 if (whitelight_instrument==
'LASCO/C1')
then
6105 lasco_rsl=5.6d0/instrument_resolution_factor
6107 else if (whitelight_instrument==
'LASCO/C2')
then
6108 lasco_rsl=11.4d0/instrument_resolution_factor
6110 else if (whitelight_instrument==
'LASCO/C3')
then
6111 lasco_rsl=56.d0/instrument_resolution_factor
6114 if (r_occultor>1.d0) r_occult=r_occultor
6115 r_occult=r_occult*const_rsun/unit_length
6116 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,ne)
6119 double precision :: nH_dummy(ixI^S)
6120 call eos%get_ne_nH(ixi^l, ixo^l, ps(igrid)%w, ne, nh_dummy)
6123 pixel=lasco_rsl*arcsec
6124 sigma0=sigma_psf*pixel
6127 r_thick=r_opt_thick*const_rsun/unit_length
6128 {
do ix^d=ixomin^d,ixomax^d\}
6129 x_sph(1:3)=ps(igrid)%x(ix^d,1:3)
6130 dx_sph(1:3)=ps(igrid)%dx(ix^d,1:3)
6132 dxl(2)=x_sph(1)*dx_sph(2)
6133 dxl(3)=x_sph(1)*dsin(x_sph(2))*dx_sph(3)
6134 ne0=ne(ix^d)*unit_numberdensity
6139 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
6141 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
6143 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
6145 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
6147 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
6149 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
6154 xsubc(1)=x_sph(1)-half*dx_sph(1)+(isubc1-half)*dx_sph(1)/nsubc1
6156 dxsubc1=dx_sph(1)/nsubc1
6160 xsubc(2)=x_sph(2)-half*dx_sph(2)+(isubc2-half)*dx_sph(2)/nsubc2
6161 dxsubc2=xsubc(1)*dx_sph(2)/nsubc2
6162 dxsubc3=xsubc(1)*dsin(xsubc(2))*dx_sph(3)/nsubc3
6163 dvolume=dxsubc1*dxsubc2*dxsubc3
6166 xsubc(3)=x_sph(3)-half*dx_sph(3)+(isubc3-half)*dx_sph(3)/nsubc3
6168 rc=dsqrt(xcent(1)**2+xcent(2)**2)
6171 if (rc>r_occult)
then
6175 tbsubc=tbsubc*dvolume*unit_length/dxi/dxi
6176 pbsubc=pbsubc*dvolume*unit_length/dxi/dxi
6177 if (tbsubc<1.d-20) emit=.false.
6182 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
6183 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
6184 ixpmin1=max(1,ixp1-3)
6185 ixpmax1=min(ixp1+3,numxi1)
6186 ixpmin2=max(1,ixp2-3)
6187 ixpmax2=min(ixp2+3,numxi2)
6188 do ixp1=ixpmin1,ixpmax1
6189 do ixp2=ixpmin2,ixpmax2
6190 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6191 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6192 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6193 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6194 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
6195 wlb(ixp1,ixp2,1)=wlb(ixp1,ixp2,1)+tbsubc*factor
6196 wlb(ixp1,ixp2,2)=wlb(ixp1,ixp2,2)+pbsubc*factor
6212 double precision,
intent(in) :: Rl
6213 double precision,
intent(inout) :: A,B,C,D
6215 double precision :: sinO,cosO,sinO2,cosO2,tmp
6220 coso=abs(dsqrt(coso2))
6221 tmp=log((1.d0+sino)/coso)
6223 b=-(1.d0-3.d0*sino2-(coso2/sino)*(1.d0+3.d0*sino2)*tmp)/8.d0
6224 c=4.d0/3.d0-coso-coso*coso2/3.d0
6225 d=(5.d0+sino2-(coso2/sino)*(5.d0-sino2)*tmp)/8.d0
6231 double precision,
intent(in) :: Rl,Rin,Ne,A,B,C,D
6232 double precision,
intent(inout) :: fluxTB,fluxPB
6234 double precision :: const,u,Bt,Br,PB,TB,sinchi2
6237 const=1.24878d-25/(1.d0-u/3.d0)
6239 bt=const*(c+u*(d-c))
6240 pb=const*sinchi2*((a+u*(b-a)))
6249 double precision,
intent(in) :: x_sph(1:3)
6250 double precision,
intent(inout) :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
6252 unitv_r(1)=dsin(x_sph(2))*dcos(x_sph(3))
6253 unitv_r(2)=dsin(x_sph(2))*dsin(x_sph(3))
6254 unitv_r(3)=dcos(x_sph(2))
6255 unitv_theta(1)=dcos(x_sph(2))*dcos(x_sph(3))
6256 unitv_theta(2)=dcos(x_sph(2))*dsin(x_sph(3))
6257 unitv_theta(3)=-dsin(x_sph(2))
6258 unitv_phi(1)=-dsin(x_sph(3))
6259 unitv_phi(2)=dcos(x_sph(3))
6264 subroutine output_data(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,datatype)
6268 integer,
intent(in) :: qunit,nXO1,nXO2,nWO
6269 double precision,
intent(in) :: dxO1(nxO1),dxO2(nxO2)
6270 double precision,
intent(in) :: xO1(nXO1),xO2(nxO2)
6271 double precision,
intent(inout) :: wO(nXO1,nXO2,nWO)
6272 character(20),
intent(in) :: datatype
6274 integer :: nPiece,nP1,nP2,nC1,nC2,nWC
6275 integer :: piece_nmax1,piece_nmax2,ix1,ix2,j,ipc,ixc1,ixc2
6276 double precision :: uniform_tol
6277 double precision,
allocatable :: xC(:,:,:,:),wC(:,:,:,:),dxC(:,:,:,:)
6284 if (abs(wo(ix1,ix2,j))<smalldouble) wo(ix1,ix2,j)=zero
6291 if (datatype==
'image_euv' .or. datatype==
'image_sxr')
then
6293 piece_nmax1=block_nx2
6294 piece_nmax2=block_nx3
6296 piece_nmax1=block_nx3
6297 piece_nmax2=block_nx1
6299 piece_nmax1=block_nx1
6300 piece_nmax2=block_nx2
6302 else if (datatype==
'spectrum_euv')
then
6305 piece_nmax2=block_nx1
6307 piece_nmax2=block_nx2
6309 piece_nmax2=block_nx3
6316 loopn1:
do j=piece_nmax1,1,-1
6317 if(mod(nxo1,j)==0)
then
6322 loopn2:
do j=piece_nmax2,1,-1
6323 if(mod(nxo2,j)==0)
then
6336 case(
'EIvtuCCmpi',
'ESvtuCCmpi',
'SIvtuCCmpi',
'WIvtuCCmpi')
6338 allocate(xc(npiece,nc1,nc2,2))
6339 allocate(dxc(npiece,nc1,nc2,2))
6340 allocate(wc(npiece,nc1,nc2,nwo))
6344 ix1=mod(ipc-1,np1)*nc1+ixc1
6345 ix2=floor(1.0*(ipc-1)/np1)*nc2+ixc2
6346 xc(ipc,ixc1,ixc2,1)=xo1(ix1)
6347 xc(ipc,ixc1,ixc2,2)=xo2(ix2)
6348 dxc(ipc,ixc1,ixc2,1)=dxo1(ix1)
6349 dxc(ipc,ixc1,ixc2,2)=dxo2(ix2)
6351 wc(ipc,ixc1,ixc2,j)=wo(ix1,ix2,j)
6358 deallocate(xc,dxc,wc)
6359 case(
'EIvtiCCmpi',
'ESvtiCCmpi',
'SIvtiCCmpi',
'WIvtiCCmpi')
6361 (maxval(abs(dxo1(:)-dxo1(1)))>uniform_tol*max(one,abs(dxo1(1))) .or. &
6362 maxval(abs(dxo2(:)-dxo2(1)))>uniform_tol*max(one,abs(dxo2(1)))))
then
6363 call mpistop(
"vti needs uniform dat-resolution image grids")
6365 call write_image_vticc(qunit,xo1,xo2,dxo1,dxo2,wo,nxo1,nxo2,nwo,nc1,nc2)
6368 call mpistop(
"Error in synthesize emission: Unknown convert_type")
6374 subroutine write_image_vticc(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,nC1,nC2)
6378 integer,
intent(in) :: qunit,nXO1,nXO2,nWO,nC1,nC2
6379 double precision,
intent(in) :: xO1(nXO1),xO2(nxO2)
6380 double precision,
intent(in) :: dxO1(nxO1),dxO2(nxO2)
6381 double precision,
intent(in) :: wO(nXO1,nXO2,nWO)
6383 double precision :: origin(1:3), spacing(1:3)
6384 integer :: wholeExtent(1:6)
6386 integer :: ixC1,ixC2
6390 character (70) :: subname,wname,vname,nameL,nameS
6391 character (len=std_len) :: filename
6392 logical :: sph_datres_no_doppler
6395 origin(1)=xo1(1)-0.5d0*dxo1(1)
6396 origin(2)=xo2(1)-0.5d0*dxo2(1)
6407 inquire(qunit,opened=fileopen)
6408 if(.not.fileopen)
then
6413 write(filename,
'(a,i4.4,a)') trim(
filename_euv),filenr,
".vti"
6415 write(filename,
'(a,i4.4,a)') trim(
filename_sxr),filenr,
".vti"
6421 open(qunit,file=filename,status=
'unknown',form=
'formatted')
6425 write(qunit,
'(a)')
'<?xml version="1.0"?>'
6426 write(qunit,
'(a)',advance=
'no')
'<VTKFile type="ImageData"'
6427 write(qunit,
'(a)')
' version="0.1" byte_order="LittleEndian">'
6428 write(qunit,
'(a,3(1pe14.6),a,6(i10),a,3(1pe14.6),a)')
' <ImageData Origin="',&
6429 origin,
'" WholeExtent="',wholeextent,
'" Spacing="',spacing,
'">'
6431 write(qunit,
'(a)')
'<FieldData>'
6432 write(qunit,
'(2a)')
'<DataArray type="Float32" Name="TIME" ',&
6433 'NumberOfTuples="1" format="ascii">'
6435 write(qunit,
'(a)')
'</DataArray>'
6436 write(qunit,
'(a)')
'</FieldData>'
6438 write(qunit,
'(a,6(i10),a)')
'<Piece Extent="',wholeextent,
'">'
6439 write(qunit,
'(a)')
'<CellData>'
6450 if (trim(
emission_model)==
'pseudo_current' .and. iw==1) vname=
'pseudo_current'
6451 if (trim(
emission_model)==
'radio_ff' .and. iw==1) vname=
'radio_brightness_temperature'
6453 if (iw==2 .and.
dat_resolution .and. (.not. sph_datres_no_doppler) .and. &
6459 ((
dat_resolution .and. ((sph_datres_no_doppler .and. iw==2) .or. &
6460 ((.not. sph_datres_no_doppler) .and. iw==3))) .or. &
6474 vname=
'absorption_fraction'
6485 if (iw==1)
write(vname,
'(a)')
'B'
6486 if (iw==2)
write(vname,
'(a)')
'pB'
6494 write(qunit,
'(a,a,a)')&
6495 '<DataArray type="Float64" Name="',trim(vname),
'" format="ascii">'
6496 write(qunit,
'(200(1pe14.6))') ((wo(ixc1,ixc2,iw),ixc1=1,nxo1),ixc2=1,nxo2)
6497 write(qunit,
'(a)')
'</DataArray>'
6499 write(qunit,
'(a)')
'</CellData>'
6500 write(qunit,
'(a)')
'</Piece>'
6502 write(qunit,
'(a)')
'</ImageData>'
6503 write(qunit,
'(a)')
'</VTKFile>'
6513 integer,
intent(in) :: qunit
6514 integer,
intent(in) :: nPiece,nC1,nC2,nWC
6515 double precision,
intent(in) :: xC(nPiece,nC1,nC2,2),dxC(nPiece,nc1,nc2,2)
6516 double precision,
intent(in) :: wC(nPiece,nC1,nC2,nWC)
6517 character(20),
intent(in) :: datatype
6520 double precision :: xP(nPiece,nC1+1,nC2+1,2)
6523 character (70) :: subname,wname,vname,nameL,nameS
6524 character (len=std_len) :: filename
6525 integer :: ixC1,ixC2,ixP,ix1,ix2,j
6526 integer :: nc,np,icel,VTK_type
6527 logical :: sph_datres_no_doppler
6538 if (ix1<np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1,1,1)-0.5d0*dxc(ixp,ix1,1,1)
6539 if (ix1==np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1-1,1,1)+0.5d0*dxc(ixp,ix1-1,1,1)
6540 if (ix2<np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2,2)-0.5d0*dxc(ixp,1,ix2,2)
6541 if (ix2==np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2-1,2)+0.5d0*dxc(ixp,1,ix2-1,2)
6546 inquire(qunit,opened=fileopen)
6547 if(.not.fileopen)
then
6551 if (datatype==
'image_euv')
then
6552 write(filename,
'(a,i4.4,a)') trim(
filename_euv),filenr,
".vtu"
6553 else if (datatype==
'image_sxr')
then
6554 write(filename,
'(a,i4.4,a)') trim(
filename_sxr),filenr,
".vtu"
6555 else if (datatype==
'image_whitelight')
then
6557 else if (datatype==
'spectrum_euv')
then
6560 open(qunit,file=filename,status=
'unknown',form=
'formatted')
6563 write(qunit,
'(a)')
'<?xml version="1.0"?>'
6564 write(qunit,
'(a)',advance=
'no')
'<VTKFile type="UnstructuredGrid"'
6565 write(qunit,
'(a)')
' version="0.1" byte_order="LittleEndian">'
6566 write(qunit,
'(a)')
'<UnstructuredGrid>'
6567 write(qunit,
'(a)')
'<FieldData>'
6568 write(qunit,
'(2a)')
'<DataArray type="Float32" Name="TIME" ',&
6569 'NumberOfTuples="1" format="ascii">'
6571 write(qunit,
'(a)')
'</DataArray>'
6572 write(qunit,
'(a)')
'</FieldData>'
6574 write(qunit,
'(a,i7,a,i7,a)') &
6575 '<Piece NumberOfPoints="',np,
'" NumberOfCells="',nc,
'">'
6576 write(qunit,
'(a)')
'<CellData>'
6578 if (datatype==
'image_euv')
then
6587 if (trim(
emission_model)==
'pseudo_current') vname=
'pseudo_current'
6588 if (trim(
emission_model)==
'radio_ff') vname=
'radio_brightness_temperature'
6591 if (j==2 .and.
dat_resolution .and. (.not. sph_datres_no_doppler) .and. &
6597 ((
dat_resolution .and. ((sph_datres_no_doppler .and. j==2) .or. &
6598 ((.not. sph_datres_no_doppler) .and. j==3))) .or. &
6612 vname=
'absorption_fraction'
6614 else if (datatype==
'image_sxr')
then
6622 else if (datatype==
'image_whitelight')
then
6623 write(vname,
'(a)')
'whitelight'
6624 else if (datatype==
'spectrum_euv')
then
6631 write(qunit,
'(a,a,a)')&
6632 '<DataArray type="Float64" Name="',trim(vname),
'" format="ascii">'
6633 write(qunit,
'(200(1pe14.6))') ((wc(ixp,ixc1,ixc2,j),ixc1=1,nc1),ixc2=1,nc2)
6634 write(qunit,
'(a)')
'</DataArray>'
6636 write(qunit,
'(a)')
'</CellData>'
6637 write(qunit,
'(a)')
'<Points>'
6638 write(qunit,
'(a)')
'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
6643 write(qunit,
'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
6645 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
6647 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
6651 write(qunit,
'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
6653 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
6655 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
6658 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
6662 write(qunit,
'(a)')
'</DataArray>'
6663 write(qunit,
'(a)')
'</Points>'
6665 write(qunit,
'(a)')
'<Cells>'
6666 write(qunit,
'(a)')
'<DataArray type="Int32" Name="connectivity" format="ascii">'
6669 write(qunit,
'(4(i7))') ix1-1+(ix2-1)*np1,ix1+(ix2-1)*np1,&
6670 ix1-1+ix2*np1,ix1+ix2*np1
6673 write(qunit,
'(a)')
'</DataArray>'
6675 write(qunit,
'(a)')
'<DataArray type="Int32" Name="offsets" format="ascii">'
6677 write(qunit,
'(i7)') icel*(2**2)
6679 write(qunit,
'(a)')
'</DataArray>'
6681 write(qunit,
'(a)')
'<DataArray type="Int32" Name="types" format="ascii">'
6685 write(qunit,
'(i2)') vtk_type
6687 write(qunit,
'(a)')
'</DataArray>'
6688 write(qunit,
'(a)')
'</Cells>'
6689 write(qunit,
'(a)')
'</Piece>'
6691 write(qunit,
'(a)')
'</UnstructuredGrid>'
6692 write(qunit,
'(a)')
'</VTKFile>'
6698 double precision,
intent(in) :: vec1(1:3),vec2(1:3)
6699 double precision,
intent(out) :: res
6701 res=vec1(1)*vec2(1)+vec1(2)*vec2(2)+vec1(3)*vec2(3)
6706 double precision,
intent(in) :: vec_in1(1:3),vec_in2(1:3)
6707 double precision,
intent(out) :: vec_out(1:3)
6709 vec_out(1)=vec_in1(2)*vec_in2(3)-vec_in1(3)*vec_in2(2)
6710 vec_out(2)=vec_in1(3)*vec_in2(1)-vec_in1(1)*vec_in2(3)
6711 vec_out(3)=vec_in1(1)*vec_in2(2)-vec_in1(2)*vec_in2(1)
6717 double precision :: LOS_psi
6718 double precision :: vec_car(1:3),vec_z(1:3),vec_temp1(1:3),vec_temp2(1:3)
6719 double precision :: vec_LOS_sph(1:3),vec_xI1_sph(1:3),vec_xI2_sph(1:3)
6723 vec_los(2)=dpi*los_theta/180.d0
6734 if (los_theta==zero)
then
6736 vec_temp1(2)=dpi/2.d0
6737 vec_temp1(3)=dpi*los_phi/180.d0
6751 los_psi=dpi*image_rotate/180.d0
6752 vec_xi1=vec_temp1*cos(los_psi)-vec_temp2*sin(los_psi)
6753 vec_xi2=vec_temp2*cos(los_psi)+vec_temp1*sin(los_psi)
6764 vec_los_sph(2:3)=vec_los_sph(2:3)*180.d0/dpi
6765 vec_xi1_sph(2:3)=vec_xi1_sph(2:3)*180.d0/dpi
6766 vec_xi2_sph(2:3)=vec_xi2_sph(2:3)*180.d0/dpi
6768 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),
']'
6769 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),
']'
6770 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),
']'
6776 double precision,
intent(in) :: vec_sph(1:3)
6777 double precision,
intent(inout) :: vec_car(1:3)
6779 vec_car(1)=vec_sph(1)*dsin(vec_sph(2))*dcos(vec_sph(3))
6780 vec_car(2)=vec_sph(1)*dsin(vec_sph(2))*dsin(vec_sph(3))
6781 vec_car(3)=vec_sph(1)*dcos(vec_sph(2))
6787 double precision,
intent(in) :: vec_car(1:3)
6788 double precision,
intent(inout) :: vec_sph(1:3)
6790 vec_sph(1)=dsqrt(vec_car(1)**2+vec_car(2)**2+vec_car(3)**2)
6791 vec_sph(2)=dacos(vec_car(3)/vec_sph(1))
6792 vec_sph(3)=atan2(vec_car(2),vec_car(1))
6798 double precision :: LOS_psi
6799 double precision :: vec_z(1:3),vec_temp1(1:3),vec_temp2(1:3)
6802 vec_los(1)=-cos(dpi*los_phi/180.d0)*sin(dpi*los_theta/180.d0)
6803 vec_los(2)=-sin(dpi*los_phi/180.d0)*sin(dpi*los_theta/180.d0)
6804 vec_los(3)=-cos(dpi*los_theta/180.d0)
6810 if (los_theta==zero)
then
6811 vec_xi1(1)=cos(dpi*los_phi/180.d0)
6812 vec_xi1(2)=sin(dpi*los_phi/180.d0)
6820 los_psi=dpi*image_rotate/180.d0
6821 vec_xi1=vec_temp1*cos(los_psi)-vec_temp2*sin(los_psi)
6822 vec_xi2=vec_temp2*cos(los_psi)+vec_temp1*sin(los_psi)
6836 double precision,
intent(in) :: x_3D_sph(1:3)
6837 double precision,
intent(inout) :: x_image(1:2)
6838 double precision :: res,res_origin
6839 double precision :: x_3D(1:3)
6850 double precision,
intent(in) :: x_3D(1:3)
6851 double precision,
intent(inout) :: x_image(1:2)
6852 double precision :: res,res_origin
6856 x_image(1)=res-res_origin
6859 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
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
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
character(len=std_len) filename_sxr
Base file name for synthetic SXR emission output.
integer spectrum_wl
wave length for spectrum
integer ixghi
Upper index of grid block arrays.
logical activate_unit_arcsec
use arcsec as length unit of images/spectra
character(len=std_len) filename_spectrum
Base file name for synthetic EUV spectrum output.
double precision global_time
The global simulation time.
logical output_absorption_fraction
output absorption fraction for thick/thin EUV synthesis when available
double precision radio_beam_fwhm
Gaussian radio beam full width at half maximum in arcsec.
integer snapshotini
Resume from the snapshot with this index.
character(len=std_len) filename_euv
Base file name for synthetic EUV emission output.
logical instrument_postprocess
Post-process dat-resolution EUV images onto the instrument pixel grid.
character(len=std_len) filename_whitelight
Base file name for synthetic white light.
character(len=std_len) convert_type
Which format to use when converting.
integer, parameter rpxmin
double precision unit_length
Physical scaling factor for length.
double precision location_slit
location of the slit
double precision time_convert_factor
Conversion factor for time unit.
integer icomm
The MPI communicator.
character(len=std_len) whitelight_instrument
white light observation instrument
integer mype
The rank of the current MPI task.
double precision radio_frequency
Observing frequency for radio free-free synthesis in Hz.
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
double precision, dimension(:), allocatable, parameter d
logical slab
Cartesian geometry or not.
double precision radio_beam_pixel_size
Output pixel size for radio beam post-processing in arcsec; <=0 uses FWHM/3.
integer snapshotnext
IO: snapshot and collapsed views output numbers/labels.
logical dat_resolution
resolution of the images
double precision r_occultor
the white light emission below it (unit=Rsun) is not visible
integer, dimension(ndim) nstretchedblocks_baselevel
(even) number of (symmetrically) stretched blocks at AMR level 1, per dimension
integer npe
The number of MPI tasks.
logical output_tau
output optical-depth map for synthetic emission when available
double precision, dimension(^nd) qstretch_baselevel
stretch factor between cells at AMR level 1, per dimension
double precision unit_velocity
Physical scaling factor for velocity.
integer radsyn_segment_batch_factor
Maximum ray segments per pixel batch, as a factor of radsyn_pixel_batch; <=0 uses memory budget....
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision los_theta
direction of the line of sight (LOS)
character(len=std_len) 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 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
subroutine get_sph_intersection_datresol_spacing(dxi)
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 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 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)
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_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)