17 double precision ::
t_aia(1:101)
20 double precision ::
f_335(1:101)
33 data t_aia / 4. , 4.05, 4.1, 4.15, 4.2, 4.25, 4.3, 4.35, &
34 4.4, 4.45, 4.5, 4.55, 4.6, 4.65, 4.7, 4.75, &
35 4.8, 4.85, 4.9, 4.95, 5. , 5.05, 5.1, 5.15, &
36 5.2, 5.25, 5.3, 5.35, 5.4, 5.45, 5.5, 5.55, &
37 5.6, 5.65, 5.7, 5.75, 5.8, 5.85, 5.9, 5.95, &
38 6. , 6.05, 6.1, 6.15, 6.2, 6.25, 6.3, 6.35, &
39 6.4, 6.45, 6.5, 6.55, 6.6, 6.65, 6.7, 6.75, &
40 6.8, 6.85, 6.9, 6.95, 7. , 7.05, 7.1, 7.15, &
41 7.2, 7.25, 7.3, 7.35, 7.4, 7.45, 7.5, 7.55, &
42 7.6, 7.65, 7.7, 7.75, 7.8, 7.85, 7.9, 7.95, &
43 8. , 8.05, 8.1, 8.15, 8.2, 8.25, 8.3, 8.35, &
44 8.4, 8.45, 8.5, 8.55, 8.6, 8.65, 8.7, 8.75, &
45 8.8, 8.85, 8.9, 8.95, 9. /
47 data f_94 / 4.25022959
d-37, 4.35880298
d-36, 3.57054296
d-35, 2.18175426
d-34, &
48 8.97592571
d-34, 2.68512961
d-33, 7.49559346
d-33, 2.11603751
d-32, &
49 5.39752853
d-32, 1.02935904
d-31, 1.33822307
d-31, 1.40884290
d-31, &
50 1.54933156
d-31, 2.07543102
d-31, 3.42026227
d-31, 6.31171444
d-31, &
51 1.16559416
d-30, 1.95360497
d-30, 2.77818735
d-30, 3.43552578
d-30, &
52 4.04061803
d-30, 4.75470982
d-30, 5.65553769
d-30, 6.70595782
d-30, &
53 7.80680354
d-30, 8.93247715
d-30, 1.02618156
d-29, 1.25979030
d-29, &
54 1.88526483
d-29, 3.62448572
d-29, 7.50553279
d-29, 1.42337571
d-28, &
55 2.37912813
d-28, 3.55232305
d-28, 4.84985757
d-28, 6.20662827
d-28, &
56 7.66193687
d-28, 9.30403645
d-28, 1.10519802
d-27, 1.25786927
d-27, &
57 1.34362634
d-27, 1.33185242
d-27, 1.22302081
d-27, 1.05677973
d-27, &
58 9.23064720
d-28, 8.78570994
d-28, 8.02397416
d-28, 5.87681142
d-28, &
59 3.82272695
d-28, 3.11492649
d-28, 3.85736090
d-28, 5.98893519
d-28, &
60 9.57553548
d-28, 1.46650267
d-27, 2.10365847
d-27, 2.79406671
d-27, &
61 3.39420087
d-27, 3.71077520
d-27, 3.57296767
d-27, 2.95114380
d-27, &
62 2.02913103
d-27, 1.13361825
d-27, 5.13405629
d-28, 2.01305089
d-28, &
63 8.15781482
d-29, 4.28366817
d-29, 3.08701543
d-29, 2.68693906
d-29, &
64 2.51764203
d-29, 2.41773103
d-29, 2.33996083
d-29, 2.26997246
d-29, &
65 2.20316143
d-29, 2.13810001
d-29, 2.07424438
d-29, 2.01149189
d-29, &
66 1.94980213
d-29, 1.88917920
d-29, 1.82963583
d-29, 1.77116920
d-29, &
67 1.71374392
d-29, 1.65740593
d-29, 1.60214447
d-29, 1.54803205
d-29, &
68 1.49510777
d-29, 1.44346818
d-29, 1.39322305
d-29, 1.34441897
d-29, &
69 1.29713709
d-29, 1.25132618
d-29, 1.20686068
d-29, 1.14226584
d-29, &
70 1.09866413
d-29, 1.05635524
d-29, 1.01532444
d-29, 9.75577134
d-30, &
71 9.37102736
d-30, 8.99873335
d-30, 8.63860172
d-30, 8.29051944
d-30, &
74 data f_131 / 3.18403601
d-37, 3.22254703
d-36, 2.61657920
d-35, &
75 1.59575286
d-34, 6.65779556
d-34, 2.07015132
d-33, &
76 6.05768615
d-33, 1.76074833
d-32, 4.52633001
d-32, &
77 8.57121883
d-32, 1.09184271
d-31, 1.10207963
d-31, &
78 1.11371658
d-31, 1.29105226
d-31, 1.80385897
d-31, &
79 3.27295431
d-31, 8.92002136
d-31, 3.15214579
d-30, &
80 9.73440787
d-30, 2.22709702
d-29, 4.01788984
d-29, &
81 6.27471832
d-29, 8.91764995
d-29, 1.18725647
d-28, &
82 1.52888040
d-28, 2.05082946
d-28, 3.47651873
d-28, &
83 8.80482184
d-28, 2.66533063
d-27, 7.05805149
d-27, &
84 1.46072515
d-26, 2.45282476
d-26, 3.55303726
d-26, &
85 4.59075911
d-26, 5.36503515
d-26, 5.68444094
d-26, &
86 5.47222296
d-26, 4.81119761
d-26, 3.85959059
d-26, &
87 2.80383406
d-26, 1.83977650
d-26, 1.11182849
d-26, &
88 6.50748885
d-27, 3.96843481
d-27, 2.61876319
d-27, &
89 1.85525324
d-27, 1.39717024
d-27, 1.11504283
d-27, &
90 9.38169611
d-28, 8.24801234
d-28, 7.43331919
d-28, &
91 6.74537063
d-28, 6.14495760
d-28, 5.70805277
d-28, &
92 5.61219786
d-28, 6.31981777
d-28, 9.19747307
d-28, &
93 1.76795732
d-27, 3.77985446
d-27, 7.43166191
d-27, &
94 1.19785603
d-26, 1.48234676
d-26, 1.36673114
d-26, &
95 9.61047146
d-27, 5.61209353
d-27, 3.04779780
d-27, &
96 1.69378976
d-27, 1.02113491
d-27, 6.82223774
d-28, &
97 5.02099099
d-28, 3.99377760
d-28, 3.36279037
d-28, &
98 2.94767378
d-28, 2.65740865
d-28, 2.44396277
d-28, &
99 2.28003967
d-28, 2.14941419
d-28, 2.04178995
d-28, &
100 1.95031045
d-28, 1.87011994
d-28, 1.79777869
d-28, &
101 1.73093957
d-28, 1.66795789
d-28, 1.60785455
d-28, &
102 1.55002399
d-28, 1.49418229
d-28, 1.44022426
d-28, &
103 1.38807103
d-28, 1.33772767
d-28, 1.28908404
d-28, &
104 1.24196208
d-28, 1.17437501
d-28, 1.12854330
d-28, &
105 1.08410498
d-28, 1.04112003
d-28, 9.99529904
d-29, &
106 9.59358806
d-29, 9.20512291
d-29, 8.83009123
d-29, &
107 8.46817043
d-29, 8.11921928
d-29 /
109 data f_171 / 2.98015581
d-42, 1.24696230
d-40, 3.37614652
d-39, 5.64103034
d-38, &
110 5.20550266
d-37, 2.77785939
d-36, 1.16283616
d-35, 6.50007689
d-35, &
111 9.96177399
d-34, 1.89586076
d-32, 2.10982799
d-31, 1.36946479
d-30, &
112 6.27396553
d-30, 2.29955134
d-29, 7.13430211
d-29, 1.91024282
d-28, &
113 4.35358848
d-28, 7.94807808
d-28, 1.07431875
d-27, 1.08399488
d-27, &
114 9.16212938
d-28, 7.34715770
d-28, 6.59246382
d-28, 9.13541375
d-28, &
115 2.05939035
d-27, 5.08206555
d-27, 1.10148083
d-26, 2.01884662
d-26, &
116 3.13578384
d-26, 4.14367719
d-26, 5.36067711
d-26, 8.74170213
d-26, &
117 1.64161233
d-25, 2.94587860
d-25, 4.76298332
d-25, 6.91765639
d-25, &
118 9.08825111
d-25, 1.08496183
d-24, 1.17440114
d-24, 1.13943939
d-24, &
119 9.71696981
d-25, 7.09593688
d-25, 4.31376399
d-25, 2.12708486
d-25, &
120 8.47429567
d-26, 3.17608104
d-26, 1.95898842
d-26, 1.98064242
d-26, &
121 1.67706555
d-26, 8.99126003
d-27, 3.29773878
d-27, 1.28896127
d-27, &
122 8.51169698
d-28, 7.53520167
d-28, 6.18268143
d-28, 4.30034650
d-28, &
123 2.78152409
d-28, 1.95437088
d-28, 1.65896278
d-28, 1.68740181
d-28, &
124 1.76054383
d-28, 1.63978419
d-28, 1.32880591
d-28, 1.00833205
d-28, &
125 7.82252806
d-29, 6.36181741
d-29, 5.34633869
d-29, 4.58013864
d-29, &
126 3.97833422
d-29, 3.49414760
d-29, 3.09790940
d-29, 2.76786227
d-29, &
127 2.48806269
d-29, 2.24823367
d-29, 2.04016653
d-29, 1.85977413
d-29, &
128 1.70367499
d-29, 1.56966125
d-29, 1.45570643
d-29, 1.35964565
d-29, &
129 1.27879263
d-29, 1.21016980
d-29, 1.15132499
d-29, 1.09959628
d-29, &
130 1.05307482
d-29, 1.01040261
d-29, 9.70657096
d-30, 9.33214234
d-30, &
131 8.97689427
d-30, 8.63761192
d-30, 8.31149879
d-30, 7.85162401
d-30, &
132 7.53828281
d-30, 7.23559452
d-30, 6.94341530
d-30, 6.66137038
d-30, &
133 6.38929156
d-30, 6.12669083
d-30, 5.87346434
d-30, 5.62943622
d-30, &
136 data f_193 / 6.40066486
d-32, 4.92737300
d-31, 2.95342934
d-30, 1.28061594
d-29, &
137 3.47747667
d-29, 5.88554792
d-29, 7.72171179
d-29, 9.75609282
d-29, &
138 1.34318963
d-28, 1.96252638
d-28, 2.70163878
d-28, 3.63192965
d-28, &
139 5.28087341
d-28, 8.37821446
d-28, 1.39089159
d-27, 2.31749718
d-27, &
140 3.77510689
d-27, 5.85198594
d-27, 8.26021568
d-27, 1.04870405
d-26, &
141 1.25209374
d-26, 1.47406787
d-26, 1.77174067
d-26, 2.24098537
d-26, &
142 3.05926105
d-26, 4.50018853
d-26, 6.84720216
d-26, 1.00595861
d-25, &
143 1.30759222
d-25, 1.36481773
d-25, 1.15943558
d-25, 1.01467304
d-25, &
144 1.04092532
d-25, 1.15071251
d-25, 1.27416033
d-25, 1.38463476
d-25, &
145 1.47882726
d-25, 1.57041238
d-25, 1.69786224
d-25, 1.94970397
d-25, &
146 2.50332918
d-25, 3.58321431
d-25, 5.18061550
d-25, 6.60405549
d-25, &
147 6.64085365
d-25, 4.83825816
d-25, 2.40545020
d-25, 8.59534098
d-26, &
148 2.90920638
d-26, 1.33204845
d-26, 9.03933926
d-27, 7.78910836
d-27, &
149 7.29342321
d-27, 7.40267022
d-27, 8.05279981
d-27, 8.13829291
d-27, &
150 6.92634262
d-27, 5.12521880
d-27, 3.59527615
d-27, 2.69617560
d-27, &
151 2.84432713
d-27, 5.06697306
d-27, 1.01281903
d-26, 1.63526978
d-26, &
152 2.06759342
d-26, 2.19482312
d-26, 2.10050611
d-26, 1.89837248
d-26, &
153 1.66347131
d-26, 1.43071097
d-26, 1.21518419
d-26, 1.02078343
d-26, &
154 8.46936184
d-27, 6.93015742
d-27, 5.56973237
d-27, 4.38951754
d-27, &
155 3.38456457
d-27, 2.55309556
d-27, 1.88904224
d-27, 1.38057546
d-27, &
156 1.00718330
d-27, 7.43581116
d-28, 5.63562931
d-28, 4.43359435
d-28, &
157 3.63923535
d-28, 3.11248143
d-28, 2.75586846
d-28, 2.50672237
d-28, &
158 2.32419348
d-28, 2.18325682
d-28, 2.06834486
d-28, 1.93497044
d-28, &
159 1.84540751
d-28, 1.76356504
d-28, 1.68741425
d-28, 1.61566157
d-28, &
160 1.54754523
d-28, 1.48249410
d-28, 1.42020176
d-28, 1.36045230
d-28, &
163 data f_211 / 4.74439912
d-42, 1.95251522
d-40, 5.19700194
d-39, 8.53120166
d-38, &
164 7.72745727
d-37, 4.04158559
d-36, 1.64853511
d-35, 8.56295439
d-35, &
165 1.17529722
d-33, 2.16867729
d-32, 2.40472264
d-31, 1.56418133
d-30, &
166 7.20032889
d-30, 2.65838271
d-29, 8.33196904
d-29, 2.26128236
d-28, &
167 5.24295811
d-28, 9.77791121
d-28, 1.35913489
d-27, 1.43957785
d-27, &
168 1.37591544
d-27, 1.49029886
d-27, 2.06183401
d-27, 3.31440622
d-27, &
169 5.42497318
d-27, 8.41100374
d-27, 1.17941366
d-26, 1.49269794
d-26, &
170 1.71506074
d-26, 1.71266353
d-26, 1.51434781
d-26, 1.36766622
d-26, &
171 1.33483562
d-26, 1.36834518
d-26, 1.45829002
d-26, 1.62575306
d-26, &
172 1.88773347
d-26, 2.22026986
d-26, 2.54930499
d-26, 2.80758138
d-26, &
173 3.06176409
d-26, 3.62799792
d-26, 5.13226109
d-26, 8.46260744
d-26, &
174 1.38486586
d-25, 1.86192535
d-25, 1.78007934
d-25, 1.16548409
d-25, &
175 5.89293257
d-26, 2.69952884
d-26, 1.24891081
d-26, 6.41273176
d-27, &
176 4.08282914
d-27, 3.26463328
d-27, 2.76230280
d-27, 2.08986882
d-27, &
177 1.37658470
d-27, 8.48489381
d-28, 5.19304217
d-28, 3.19312514
d-28, &
178 2.02968197
d-28, 1.50171666
d-28, 1.39164218
d-28, 1.42448821
d-28, &
179 1.41714519
d-28, 1.33341059
d-28, 1.20759270
d-28, 1.07259692
d-28, &
180 9.44895400
d-29, 8.29030041
d-29, 7.25440631
d-29, 6.33479483
d-29, &
181 5.51563757
d-29, 4.79002469
d-29, 4.14990482
d-29, 3.59384972
d-29, &
182 3.12010860
d-29, 2.72624742
d-29, 2.40734791
d-29, 2.15543565
d-29, &
183 1.95921688
d-29, 1.80682882
d-29, 1.68695662
d-29, 1.59020936
d-29, &
184 1.50940886
d-29, 1.43956179
d-29, 1.37731622
d-29, 1.32049043
d-29, &
185 1.26771875
d-29, 1.21803879
d-29, 1.17074716
d-29, 1.10507836
d-29, &
186 1.06022834
d-29, 1.01703080
d-29, 9.75436986
d-30, 9.35349257
d-30, &
187 8.96744546
d-30, 8.59527489
d-30, 8.23678940
d-30, 7.89144480
d-30, &
190 data f_304 / 3.62695850
d-32, 2.79969087
d-31, 1.68340584
d-30, 7.32681440
d-30, &
191 1.99967770
d-29, 3.41296785
d-29, 4.55409104
d-29, 5.94994635
d-29, &
192 8.59864963
d-29, 1.39787633
d-28, 3.17701965
d-28, 1.14474920
d-27, &
193 4.44845958
d-27, 1.54785841
d-26, 4.70265345
d-26, 1.24524365
d-25, &
194 2.81535352
d-25, 5.10093666
d-25, 6.83545307
d-25, 6.82110329
d-25, &
195 5.66886188
d-25, 4.36205513
d-25, 3.29265688
d-25, 2.49802368
d-25, &
196 1.92527113
d-25, 1.51058572
d-25, 1.20596047
d-25, 9.76884267
d-26, &
197 7.89979266
d-26, 6.18224289
d-26, 4.67298332
d-26, 3.57934505
d-26, &
198 2.84535785
d-26, 2.32853022
d-26, 1.95228514
d-26, 1.67880071
d-26, &
199 1.47608785
d-26, 1.32199691
d-26, 1.20070960
d-26, 1.09378177
d-26, &
200 1.00031730
d-26, 9.62434001
d-27, 1.05063954
d-26, 1.27267143
d-26, &
201 1.45923057
d-26, 1.36746707
d-26, 1.03466970
d-26, 6.97647829
d-27, &
202 4.63141039
d-27, 3.19031994
d-27, 2.33373613
d-27, 1.81589079
d-27, &
203 1.48446917
d-27, 1.26611478
d-27, 1.12617468
d-27, 1.03625148
d-27, &
204 9.61400595
d-28, 8.79016231
d-28, 7.82612130
d-28, 6.73762960
d-28, &
205 5.59717956
d-28, 4.53010243
d-28, 3.65712196
d-28, 3.00958686
d-28, &
206 2.54011502
d-28, 2.18102277
d-28, 1.88736437
d-28, 1.63817539
d-28, &
207 1.42283147
d-28, 1.23631916
d-28, 1.07526003
d-28, 9.36797928
d-29, &
208 8.18565660
d-29, 7.18152734
d-29, 6.32523238
d-29, 5.59513985
d-29, &
209 4.96614048
d-29, 4.42518826
d-29, 3.95487628
d-29, 3.54690294
d-29, &
210 3.18953930
d-29, 2.87720933
d-29, 2.60186750
d-29, 2.36011522
d-29, &
211 2.14717806
d-29, 1.95905217
d-29, 1.79287981
d-29, 1.64562262
d-29, &
212 1.51489425
d-29, 1.39876064
d-29, 1.29496850
d-29, 1.18665438
d-29, &
213 1.10240474
d-29, 1.02643099
d-29, 9.57780996
d-30, 8.95465151
d-30, &
214 8.38950190
d-30, 7.87283711
d-30, 7.40136507
d-30, 6.96804279
d-30, &
217 data f_335 / 2.46882661
d-32, 1.89476632
d-31, 1.13216502
d-30, 4.89532008
d-30, &
218 1.32745970
d-29, 2.25390335
d-29, 3.00511672
d-29, 3.96035934
d-29, &
219 5.77977656
d-29, 8.58600736
d-29, 1.14083000
d-28, 1.48644411
d-28, &
220 2.15788823
d-28, 3.51628877
d-28, 6.12200698
d-28, 1.08184987
d-27, &
221 1.85590697
d-27, 2.91679107
d-27, 3.94405396
d-27, 4.63610680
d-27, &
222 5.13824456
d-27, 5.66602209
d-27, 6.30009232
d-27, 7.03422868
d-27, &
223 7.77973918
d-27, 8.32371831
d-27, 8.56724316
d-27, 8.62601374
d-27, &
224 8.13308844
d-27, 6.53188216
d-27, 4.55197029
d-27, 3.57590087
d-27, &
225 3.59571707
d-27, 4.03502770
d-27, 4.54366411
d-27, 4.96914990
d-27, &
226 5.24601170
d-27, 5.39979250
d-27, 5.43023669
d-27, 5.26235042
d-27, &
227 4.91585495
d-27, 4.52628362
d-27, 4.13385020
d-27, 3.67538967
d-27, &
228 3.39939742
d-27, 3.81284533
d-27, 5.02332701
d-27, 6.19438602
d-27, &
229 6.49613071
d-27, 6.04010475
d-27, 5.24664275
d-27, 4.37225997
d-27, &
230 3.52957182
d-27, 2.76212276
d-27, 2.08473158
d-27, 1.50850518
d-27, &
231 1.04602472
d-27, 7.13091243
d-28, 5.34289645
d-28, 5.21079581
d-28, &
232 6.22246365
d-28, 6.99555864
d-28, 6.29665489
d-28, 4.45077026
d-28, &
233 2.67046793
d-28, 1.52774686
d-28, 9.18061770
d-29, 6.09116074
d-29, &
234 4.48562572
d-29, 3.59463696
d-29, 3.05820218
d-29, 2.70766652
d-29, &
235 2.46144034
d-29, 2.27758450
d-29, 2.13331183
d-29, 2.01537836
d-29, &
236 1.91566180
d-29, 1.82893912
d-29, 1.75167748
d-29, 1.68136168
d-29, &
237 1.61615595
d-29, 1.55481846
d-29, 1.49643236
d-29, 1.44046656
d-29, &
238 1.38657085
d-29, 1.33459068
d-29, 1.28447380
d-29, 1.23615682
d-29, &
239 1.18963296
d-29, 1.14478976
d-29, 1.10146637
d-29, 1.04039479
d-29, &
240 9.98611410
d-30, 9.58205147
d-30, 9.19202009
d-30, 8.81551313
d-30, &
241 8.45252127
d-30, 8.10224764
d-30, 7.76469090
d-30, 7.43954323
d-30, &
247 data t_iris / 4. , 4.1 , 4.2 , 4.3 , 4.40000001, &
248 4.50000001, 4.60000001, 4.70000001, 4.80000001, 4.90000001, &
249 5.00000001, 5.10000002, 5.20000002, 5.30000002, 5.40000002, &
250 5.50000002, 5.60000002, 5.70000003, 5.80000003, 5.90000003, &
251 6.00000003, 6.10000003, 6.20000003, 6.30000003, 6.40000004, &
252 6.50000004, 6.60000004, 6.70000004, 6.80000004, 6.90000004, &
253 7.00000004, 7.10000005, 7.20000005, 7.30000005, 7.40000005, &
254 7.50000005, 7.60000005, 7.70000006, 7.80000006, 7.90000006, &
257 data f_1354 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
258 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
259 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, 1.09503647
d-39, &
263 5.47214550
d-36, 2.42433983
d-33, 2.75295034
d-31, 1.21929718
d-29, &
264 2.48392125
d-28, 2.33268145
d-27, 8.68623633
d-27, 1.00166284
d-26, &
265 3.63126633
d-27, 7.45174807
d-28, 1.38224064
d-28, 2.69270994
d-29, &
266 5.53314977
d-30, 1.15313092
d-30, 2.34195788
d-31, 4.48242942
d-32, &
272 data t_eis1 / 1.99526231
d+05, 2.23872114
d+05, 2.51188643
d+05, 2.81838293
d+05, &
273 3.16227766
d+05, 3.54813389
d+05, 3.98107171
d+05, 4.46683592
d+05, &
274 5.01187234
d+05, 5.62341325
d+05, 6.30957344
d+05, 7.07945784
d+05, &
275 7.94328235
d+05, 8.91250938
d+05, 1.00000000
d+06, 1.12201845
d+06, &
276 1.25892541
d+06, 1.41253754
d+06, 1.58489319
d+06, 1.77827941
d+06, &
277 1.99526231
d+06, 2.23872114
d+06, 2.51188643
d+06, 2.81838293
d+06, &
278 3.16227766
d+06, 3.54813389
d+06, 3.98107171
d+06, 4.46683592
d+06, &
279 5.01187234
d+06, 5.62341325
d+06, 6.30957344
d+06, 7.07945784
d+06, &
280 7.94328235
d+06, 8.91250938
d+06, 1.00000000
d+07, 1.12201845
d+07, &
281 1.25892541
d+07, 1.41253754
d+07, 1.58489319
d+07, 1.77827941
d+07, &
282 1.99526231
d+07, 2.23872114
d+07, 2.51188643
d+07, 2.81838293
d+07, &
283 3.16227766
d+07, 3.54813389
d+07, 3.98107171
d+07, 4.46683592
d+07, &
284 5.01187234
d+07, 5.62341325
d+07, 6.30957344
d+07, 7.07945784
d+07, &
285 7.94328235
d+07, 8.91250938
d+07, 1.00000000
d+08, 1.12201845
d+08, &
286 1.25892541
d+08, 1.41253754
d+08, 1.58489319
d+08, 1.77827941
d+08 /
288 data t_eis2 / 1.99526231
d+06, 2.23872114
d+06, 2.51188643
d+06, 2.81838293
d+06, &
289 3.16227766
d+06, 3.54813389
d+06, 3.98107171
d+06, 4.46683592
d+06, &
290 5.01187234
d+06, 5.62341325
d+06, 6.30957344
d+06, 7.07945784
d+06, &
291 7.94328235
d+06, 8.91250938
d+06, 1.00000000
d+07, 1.12201845
d+07, &
292 1.25892541
d+07, 1.41253754
d+07, 1.58489319
d+07, 1.77827941
d+07, &
293 1.99526231
d+07, 2.23872114
d+07, 2.51188643
d+07, 2.81838293
d+07, &
294 3.16227766
d+07, 3.54813389
d+07, 3.98107171
d+07, 4.46683592
d+07, &
295 5.01187234
d+07, 5.62341325
d+07, 6.30957344
d+07, 7.07945784
d+07, &
296 7.94328235
d+07, 8.91250938
d+07, 1.00000000
d+08, 1.12201845
d+08, &
297 1.25892541
d+08, 1.41253754
d+08, 1.58489319
d+08, 1.77827941
d+08, &
298 1.99526231
d+08, 2.23872114
d+08, 2.51188643
d+08, 2.81838293
d+08, &
299 3.16227766
d+08, 3.54813389
d+08, 3.98107171
d+08, 4.46683592
d+08, &
300 5.01187234
d+08, 5.62341325
d+08, 6.30957344
d+08, 7.07945784
d+08, &
301 7.94328235
d+08, 8.91250938
d+08, 1.00000000
d+09, 1.12201845
d+09, &
302 1.25892541
d+09, 1.41253754
d+09, 1.58489319
d+09, 1.77827941
d+09 /
304 data f_263 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
305 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
306 0.00000000
d+00, 4.46454917
d-45, 3.26774829
d-42, 1.25292566
d-39, &
307 2.66922338
d-37, 3.28497742
d-35, 2.38677554
d-33, 1.03937729
d-31, &
308 2.75075687
d-30, 4.47961733
d-29, 4.46729177
d-28, 2.64862689
d-27, &
309 8.90863800
d-27, 1.72437548
d-26, 2.22217752
d-26, 2.27999477
d-26, &
310 2.08264363
d-26, 1.78226687
d-26, 1.45821699
d-26, 1.14675379
d-26, &
311 8.63082492
d-27, 6.15925429
d-27, 4.11252514
d-27, 2.51530564
d-27, &
312 1.37090986
d-27, 6.42443134
d-28, 2.48392636
d-28, 7.59187874
d-29, &
313 1.77852938
d-29, 3.23945221
d-30, 4.90533903
d-31, 6.75458158
d-32, &
314 9.06878868
d-33, 1.23927474
d-33, 1.75769395
d-34, 2.60710914
d-35, &
315 4.04318030
d-36, 6.53500581
d-37, 1.09365022
d-37, 1.88383322
d-38, &
316 3.31425233
d-39, 5.90964084
d-40, 1.06147549
d-40, 1.90706170
d-41, &
317 3.41331584
d-42, 6.07310718
d-43, 1.07364738
d-43, 1.89085498
d-44, &
318 3.32598922
d-45, 5.87125640
d-46, 0.00000000
d+00, 0.00000000
d+00 /
320 data f_264 / 0.00000000
d+00, 2.81670057
d-46, 1.28007268
d-43, 2.54586603
d-41, &
321 2.67887256
d-39, 1.68413285
d-37, 6.85702304
d-36, 1.91797284
d-34, &
322 3.84675839
d-33, 5.69939170
d-32, 6.36224608
d-31, 5.39176489
d-30, &
323 3.45478458
d-29, 1.64848693
d-28, 5.71476364
d-28, 1.39909997
d-27, &
324 2.37743056
d-27, 2.86712530
d-27, 2.65206348
d-27, 2.07175767
d-27, &
325 1.47866767
d-27, 1.01087374
d-27, 6.79605811
d-28, 4.54746770
d-28, &
326 3.04351751
d-28, 2.03639149
d-28, 1.35940991
d-28, 9.01451939
d-29, &
327 5.91289972
d-29, 3.81821178
d-29, 2.41434696
d-29, 1.48871220
d-29, &
328 8.93362094
d-30, 5.21097445
d-30, 2.95964719
d-30, 1.64278748
d-30, &
329 8.95571660
d-31, 4.82096011
d-31, 2.57390991
d-31, 1.36821781
d-31, &
330 7.27136350
d-32, 3.87019426
d-32, 2.06883430
d-32, 1.11228884
d-32, &
331 6.01883313
d-33, 3.27790676
d-33, 1.79805012
d-33, 9.93085346
d-34, &
332 5.52139556
d-34, 3.08881387
d-34, 1.73890315
d-34, 9.84434964
d-35, &
333 5.60603378
d-35, 3.20626492
d-35, 1.84111068
d-35, 0.00000000
d+00, &
334 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00 /
336 data f_192 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 4.35772105
d-44, &
337 1.26162319
d-41, 1.97471205
d-39, 1.83409019
d-37, 1.08206288
d-35, &
338 4.27914363
d-34, 1.17943846
d-32, 2.32565755
d-31, 3.33087991
d-30, &
339 3.47013260
d-29, 2.60375866
d-28, 1.37737127
d-27, 5.01053913
d-27, &
340 1.23479810
d-26, 2.11310542
d-26, 2.71831513
d-26, 2.89851163
d-26, &
341 2.77312376
d-26, 2.50025229
d-26, 2.18323661
d-26, 1.86980322
d-26, &
342 1.58035034
d-26, 1.31985651
d-26, 1.08733133
d-26, 8.81804906
d-27, &
343 7.00417973
d-27, 5.43356567
d-27, 4.09857884
d-27, 2.99651764
d-27, &
344 2.11902962
d-27, 1.45014127
d-27, 9.62291023
d-28, 6.21548647
d-28, &
345 3.92807578
d-28, 2.44230375
d-28, 1.50167782
d-28, 9.17611405
d-29, &
346 5.58707641
d-29, 3.40570915
d-29, 2.08030862
d-29, 1.27588676
d-29, &
347 7.86535588
d-30, 4.87646151
d-30, 3.03888897
d-30, 1.90578649
d-30, &
348 1.20195947
d-30, 7.61955060
d-31, 4.85602199
d-31, 3.11049969
d-31, &
349 2.00087065
d-31, 1.29223740
d-31, 8.37422008
d-32, 0.00000000
d+00, &
350 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00 /
352 data f_255 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 1.76014287
d-44, &
353 5.07057938
d-42, 7.90473970
d-40, 7.31852999
d-38, 4.30709255
d-36, &
354 1.70009061
d-34, 4.67925160
d-33, 9.21703546
d-32, 1.31918676
d-30, &
355 1.37393161
d-29, 1.03102379
d-28, 5.45694018
d-28, 1.98699648
d-27, &
356 4.90346776
d-27, 8.40524725
d-27, 1.08321456
d-26, 1.15714525
d-26, &
357 1.10905152
d-26, 1.00155023
d-26, 8.75799161
d-27, 7.50935839
d-27, &
358 6.35253533
d-27, 5.30919268
d-27, 4.37669455
d-27, 3.55185164
d-27, &
359 2.82347055
d-27, 2.19257595
d-27, 1.65589541
d-27, 1.21224987
d-27, &
360 8.58395132
d-28, 5.88163935
d-28, 3.90721447
d-28, 2.52593407
d-28, &
361 1.59739995
d-28, 9.93802874
d-29, 6.11343388
d-29, 3.73711135
d-29, &
362 2.27618743
d-29, 1.38793199
d-29, 8.48060787
d-30, 5.20305940
d-30, &
363 3.20867365
d-30, 1.99011277
d-30, 1.24064551
d-30, 7.78310544
d-31, &
364 4.91013681
d-31, 3.11338381
d-31, 1.98451675
d-31, 1.27135460
d-31, &
365 8.17917486
d-32, 5.28280497
d-32, 3.42357159
d-32, 0.00000000
d+00, &
366 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00 /
371 integer,
intent(in) :: ixI^L, ixO^L
372 double precision,
intent(in) :: w(ixI^S,nw)
373 double precision,
intent(in) :: x(ixI^S,1:ndim)
374 double precision,
intent(out):: res(ixI^S)
381 procedure(
get_subr1),
pointer,
nopass :: get_rho => null()
382 procedure(
get_subr1),
pointer,
nopass :: get_pthermal => null()
383 procedure(
get_subr1),
pointer,
nopass :: get_var_rfactor => null()
387 double precision :: rfactor = 1d0
394 subroutine get_line_info(wl,ion,mass,logTe,line_center,spatial_px,spectral_px,sigma_PSF,width_slit)
406 integer,
intent(in) :: wl
407 integer,
intent(out) :: mass
408 character(len=30),
intent(out) :: ion
409 double precision,
intent(out) :: logTe,line_center,spatial_px,spectral_px
410 double precision,
intent(out) :: sigma_PSF,width_slit
489 line_center=262.976d0
498 line_center=263.765d0
507 line_center=192.028d0
516 line_center=255.113d0
522 call mpistop(
"No information about this line")
526 subroutine get_euv(wl,ixI^L,ixO^L,w,x,fl,flux)
533 integer,
intent(in) :: wl
534 integer,
intent(in) :: ixI^L, ixO^L
535 double precision,
intent(in) :: x(ixI^S,1:ndim)
536 double precision,
intent(in) :: w(ixI^S,1:nw)
538 double precision,
intent(out) :: flux(ixI^S)
541 double precision,
allocatable :: t_table(:),f_table(:)
542 integer :: ix^D,iTt,i
543 double precision :: pth(ixI^S),Te(ixI^S),Ne(ixI^S)
544 double precision :: logT,logGT,GT
550 allocate(t_table(1:n_table))
551 allocate(f_table(1:n_table))
556 allocate(t_table(1:n_table))
557 allocate(f_table(1:n_table))
562 allocate(t_table(1:n_table))
563 allocate(f_table(1:n_table))
568 allocate(t_table(1:n_table))
569 allocate(f_table(1:n_table))
574 allocate(t_table(1:n_table))
575 allocate(f_table(1:n_table))
580 allocate(t_table(1:n_table))
581 allocate(f_table(1:n_table))
586 allocate(t_table(1:n_table))
587 allocate(f_table(1:n_table))
592 allocate(t_table(1:n_table))
593 allocate(f_table(1:n_table))
598 allocate(t_table(1:n_table))
599 allocate(f_table(1:n_table))
604 allocate(t_table(1:n_table))
605 allocate(f_table(1:n_table))
610 allocate(t_table(1:n_table))
611 allocate(f_table(1:n_table))
616 allocate(t_table(1:n_table))
617 allocate(f_table(1:n_table))
623 call mpistop(
"Unknown wavelength")
625 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
626 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
627 if(
associated(fl%get_var_Rfactor))
then
628 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
635 flux(ixo^s)=ne(ixo^s)**2
638 flux(ixo^s)=ne(ixo^s)**2
642 case(94,131,171,193,211,304,335,1354)
645 if(f_table(i) .gt. 1.d-99)
then
646 f_table(i)=log10(f_table(i))
652 {
do ix^db=ixomin^db,ixomax^db\}
654 if (logt>=t_table(1) .and. logt<=t_table(n_table))
then
656 if (logt>=t_table(itt) .and. logt<t_table(itt+1))
then
657 loggt=f_table(itt)*(logt-t_table(itt+1))/(t_table(itt)-t_table(itt+1))+&
658 f_table(itt+1)*(logt-t_table(itt))/(t_table(itt+1)-t_table(itt))
661 flux(ix^d)=flux(ix^d)*(10**(loggt))
662 if(flux(ix^d)<smalldouble) flux(ix^d)=0.d0
669 {
do ix^db=ixomin^db,ixomax^db\}
670 if (te(ix^d)>=t_table(1) .and. te(ix^d)<=t_table(n_table))
then
672 if (te(ix^d)>=t_table(itt) .and. te(ix^d)<t_table(itt+1))
then
673 gt=f_table(itt)*(te(ix^d)-t_table(itt+1))/(t_table(itt)-t_table(itt+1))+&
674 f_table(itt+1)*(te(ix^d)-t_table(itt))/(t_table(itt+1)-t_table(itt))
677 flux(ix^d)=flux(ix^d)*gt
678 if(flux(ix^d)<smalldouble) flux(ix^d)=0.d0
685 deallocate(t_table,f_table)
688 subroutine get_sxr(ixI^L,ixO^L,w,x,fl,flux,El,Eu)
695 integer,
intent(in) :: ixI^L,ixO^L
696 integer,
intent(in) :: El,Eu
697 double precision,
intent(in) :: x(ixI^S,1:ndim)
698 double precision,
intent(in) :: w(ixI^S,nw)
700 double precision,
intent(out) :: flux(ixI^S)
702 integer :: ix^D,ixO^D
704 double precision :: I0,kb,keV,dE,Ei
705 double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
706 double precision :: Ne(ixI^S),gff(ixI^S),fi(ixI^S)
707 double precision :: EM(ixI^S)
713 nume=floor((eu-el)/de)
714 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
715 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
717 if(
associated(fl%get_var_Rfactor))
then
718 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
725 em(ixo^s)=(ne(ixo^s))**2*1.d6
728 em(ixo^s)=(ne(ixo^s))**2
730 kbt(ixo^s)=kb*te(ixo^s)/kev
735 {
do ix^db=ixomin^db,ixomax^db\}
736 if (kbt(ix^d)>0.01*ei)
then
737 if(kbt(ix^d)<ei) gff(ix^d)=(kbt(ix^d)/ei)**0.4
738 fi(ix^d)=(em(ix^d)*gff(ix^d))*dexp(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
743 flux(ixo^s)=flux(ixo^s)+fi(ixo^s)*de
745 flux(ixo^s)=flux(ixo^s)*i0
752 double precision,
intent(in) :: xbox^L
754 double precision,
intent(out) :: eflux
756 double precision :: dxb^D,xb^L
757 integer :: iigrid,igrid,j
758 integer :: ixO^L,ixI^L,ix^D
759 double precision :: eflux_grid,eflux_pe
766 do iigrid=1,igridstail; igrid=igrids(iigrid);
770 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)
771 eflux_pe=eflux_pe+eflux_grid
773 call mpi_allreduce(eflux_pe,eflux,1,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
779 integer,
intent(in) :: ixI^L,ixO^L
780 double precision,
intent(in) :: x(ixI^S,1:ndim),dV(ixI^S)
781 double precision,
intent(in) :: w(ixI^S,nw)
782 double precision,
intent(in) :: xbox^L,xb^L
784 double precision,
intent(out) :: eflux_grid
786 integer :: ix^D,ixO^D,ixb^L
787 integer :: iE,numE,j,inbox
788 double precision :: I0,kb,keV,dE,Ei,El,Eu,A_cgs
789 double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
790 double precision :: Ne(ixI^S),EM(ixI^S)
791 double precision :: gff,fi,erg_SI
795 {
if (xbmin^d<xboxmax^d .and. xbmax^d>xboxmin^d) inbox=inbox+1\}
797 if (inbox==ndim)
then
799 ^d&ixbmin^d=ixomin^d;
800 ^d&ixbmax^d=ixomax^d;
801 {
if (xbmax^d>xboxmax^d) ixbmax^d=ixomax^d-ceiling((xbmax^d-xboxmax^d)/
dxlevel(^d))\}
802 {
if (xbmin^d<xboxmin^d) ixbmin^d=ceiling((xboxmin^d-xbmin^d)/
dxlevel(^d))+ixomin^d\}
809 el=const_h*const_c/(8.d0*a_cgs)/kev
810 eu=const_h*const_c/(1.d0*a_cgs)/kev
812 nume=floor((eu-el)/de)
813 call fl%get_pthermal(w,x,ixi^l,ixb^l,pth)
814 call fl%get_rho(w,x,ixi^l,ixb^l,ne)
815 if(
associated(fl%get_var_Rfactor))
then
816 call fl%get_var_Rfactor(w,x,ixi^l,ixb^l,te)
823 em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*(
unit_length*1.d2)**3
826 em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*
unit_length**3
828 kbt(ixb^s)=kb*te(ixb^s)/kev
833 {
do ix^db=ixbmin^db,ixbmax^db\}
834 if (kbt(ix^d)>1.d-2*ei)
then
835 if(kbt(ix^d)<ei)
then
836 gff=(kbt(ix^d)/ei)**0.4
840 fi=(em(ix^d)*gff)*dexp(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
841 eflux_grid=eflux_grid+fi*de*ei
845 eflux_grid=eflux_grid*kev*erg_si
854 integer,
intent(in) :: qunit
856 character(20) :: datatype
859 character (30) :: ion
860 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
861 double precision :: xslit,arcsec
865 if (
mype==0) print *,
'###################################################'
868 if (
mype==0) print *,
'Systhesizing EUV spectrum (observed by IRIS).'
869 case (263,264,192,255)
870 if (
mype==0) print *,
'Systhesizing EUV spectrum (observed by Hinode/EIS).'
872 call mpistop(
'Wrong wavelength!')
876 if (
mype==0)
write(*,
'(a,f8.3,a)')
' Wavelength: ',linecent,
' Angstrom'
877 if (
mype==0) print *,
'Unit of EUV flux: DN s^-1 pixel^-1'
878 if (
mype==0)
write(*,
'(a,f5.3,a,f5.1,a)')
' Pixel: ',wlrsl,
' Angstrom x ',spacersl*725.0,
' km'
881 call mpistop(
'Wrong spectrum window!')
884 datatype=
'spectrum_euv'
887 if (
mype==0) print *,
'Unit of wavelength: Angstrom (0.1 nm) '
888 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
894 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Width of slit: ',wslit*725.0,
' km'
897 if (
mype==0) print *,
'Unit of wavelength: Angstrom (0.1 nm) '
898 if (
mype==0) print *,
'Unit of length: arcsec (~725 km)'
899 if (
mype==0) print *,
'Direction of the slit: parallel to xI2 vector'
900 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Location of slit: xI1 = ',
location_slit,
' arcsec'
901 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Width of slit: ',wslit,
' arcsec'
904 call mpistop(
'Wrong resolution for resolution_spectrum!')
907 if (
mype==0) print *,
'###################################################'
913 integer,
intent(in) :: qunit
914 character(20),
intent(in) :: datatype
917 integer :: numWL,numXS,iwL,ixS,numWI,ix^D
918 double precision :: dwLg,dxSg,xSmin,xSmax,xScent,wLmin,wLmax
919 double precision,
allocatable :: wL(:),xS(:),dwL(:),dxS(:)
920 double precision,
allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
921 double precision :: vec_cor(1:3),xI_cor(1:2)
922 double precision :: res,r_loc,r_max
925 character (30) :: ion
926 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
927 double precision :: unitv,arcsec,RHESSI_rsl,pixel
928 integer :: iigrid,igrid,i,j,numS
929 double precision :: xLmin,xLmax,xslit
935 if (ix1==1) vec_cor(1)=xprobmin1
936 if (ix1==2) vec_cor(1)=xprobmax1
938 if (ix2==1) vec_cor(2)=xprobmin2
939 if (ix2==2) vec_cor(2)=xprobmax2
941 if (ix3==1) vec_cor(3)=xprobmin3
942 if (ix3==2) vec_cor(3)=xprobmax3
945 r_loc=r_loc+(vec_cor(2)-
x_origin(2))**2
946 r_loc=r_loc+(vec_cor(3)-
x_origin(3))**2
948 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
951 r_max=max(r_max,r_loc)
955 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
959 xsmin=min(xsmin,xi_cor(2))
960 xsmax=max(xsmax,xi_cor(2))
970 xscent=(xsmin+xsmax)/2.d0
980 numxs=ceiling((xsmax-xscent)/dxsg)
981 xsmin=xscent-numxs*dxsg
982 xsmax=xscent+numxs*dxsg
988 allocate(wl(numwl),dwl(numwl),xs(numxs),dxs(numxs))
990 allocate(wi(numwl,numxs,numwi),spectra(numwl,numxs),spectra_rc(numwl,numxs))
992 wl(iwl)=wlmin+iwl*dwlg-half*dwlg
996 xs(ixs)=xsmin+dxsg*(ixs-half)
1002 do iigrid=1,igridstail; igrid=igrids(iigrid);
1004 if (ix1==1) vec_cor(1)=
rnode(rpxmin1_,igrid)
1005 if (ix1==2) vec_cor(1)=
rnode(rpxmax1_,igrid)
1007 if (ix2==1) vec_cor(2)=
rnode(rpxmin2_,igrid)
1008 if (ix2==2) vec_cor(2)=
rnode(rpxmax2_,igrid)
1010 if (ix3==1) vec_cor(3)=
rnode(rpxmin3_,igrid)
1011 if (ix3==2) vec_cor(3)=
rnode(rpxmax3_,igrid)
1013 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1017 xlmin=min(xlmin,xi_cor(1))
1018 xlmax=max(xlmax,xi_cor(1))
1025 if (xslit>=xlmin-wslit*arcsec .and. xslit<=xlmax+wslit*arcsec)
then
1031 call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
1035 if (spectra_rc(iwl,ixs)>smalldouble)
then
1036 wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
1046 call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
1048 deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
1054 integer,
intent(in) :: igrid,numWL,numXS
1055 double precision,
intent(in) :: wL(numWL),xS(numXS)
1056 double precision,
intent(in) :: dwLg,dxSg
1057 double precision,
intent(inout) :: spectra(numWL,numXS)
1060 integer :: ixO^L,ixI^L,ix^D,ixOnew,j
1061 double precision,
allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
1062 double precision :: wlc,wlwd,res,dst_slit,xslit,arcsec
1063 double precision :: vloc(1:3),xloc(1:3),dxloc(1:3),xIloc(1:2),dxIloc(1:2)
1064 integer :: nSubC^D,iSubC^D,iwL,ixS,ixSmin,ixSmax,iwLmin,iwLmax,nwL
1065 double precision :: slit_width,dxSubC^D,xerf^L,fluxSubC
1066 double precision :: xSubC(1:3),xCent(1:2)
1069 double precision :: logTe,lineCent
1070 character (30) :: ion
1071 double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
1072 double precision :: sigma_wl,sigma_xs,factor
1083 ^d&ixomin^d=ixmlo^d\
1084 ^d&ixomax^d=ixmhi^d\
1085 ^d&iximin^d=
ixglo^d\
1086 ^d&iximax^d=
ixghi^d\
1087 allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
1090 call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
1091 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
1092 if(
associated(fl%get_var_Rfactor))
then
1093 call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
1094 te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
1096 te(ixo^s)=pth(ixo^s)/(fl%Rfactor*rho(ixo^s))
1098 {
do ix^d=ixomin^d,ixomax^d\}
1100 vloc(j)=ps(igrid)%w(ix^d,iw_mom(j))/rho(ix^d)
1108 slit_width=wslit*arcsec
1109 sigma_wl=sigma_psf*dwlg
1110 sigma_xs=sigma_psf*dxsg
1111 {
do ix^d=ixomin^d,ixomax^d\}
1112 xloc(1:3)=ps(igrid)%x(ix^d,1:3)
1113 dxloc(1:3)=ps(igrid)%dx(ix^d,1:3)
1117 if (xiloc(1)>=xslit-half*(slit_width+dxiloc(1)) .and. &
1118 xiloc(1)<=xslit+half*(slit_width+dxiloc(1)))
then
1120 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi1(^d))/(slit_width/16.d0)));
1121 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi2(^d))/(dxsg/4.d0)));
1122 ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
1125 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*
unit_length*1.d2/dxsg/dxsg
1128 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*
unit_length/dxsg/dxsg
1132 wlwd=wlwd*linecent/const_c
1134 {
do isubc^d=1,nsubc^d\}
1135 ^d&xsubc(^d)=xloc(^d)-half*dxloc(^d)+(isubc^d-half)*dxsubc^d;
1137 dst_slit=abs(xcent(1)-xslit)
1138 if (dst_slit<=half*slit_width)
then
1139 ixs=floor((xcent(2)-(xs(1)-half*dxsg))/dxsg)+1
1141 ixsmax=min(ixs+3,numxs)
1142 iwl=floor((wlc-(wl(1)-half*dwlg))/dwlg)+1
1143 nwl=3*ceiling(wlwd/dwlg+1)
1144 iwlmin=max(1,iwl-nwl)
1145 iwlmax=min(iwl+nwl,numwl)
1147 do iwl=iwlmin,iwlmax
1148 do ixs=ixsmin,ixsmax
1149 xerfmin1=(wl(iwl)-half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
1150 xerfmax1=(wl(iwl)+half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
1151 xerfmin2=(xs(ixs)-half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
1152 xerfmax2=(xs(ixs)+half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
1153 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
1154 spectra(iwl,ixs)=spectra(iwl,ixs)+fluxsubc*factor
1163 deallocate(flux,v,pth,te)
1168 integer,
intent(in) :: qunit
1169 character(20),
intent(in) :: datatype
1172 integer :: numWL,numXS,iwL,ixS,numWI,numS
1173 double precision :: dwLg,xSmin,xSmax,wLmin,wLmax
1174 double precision,
allocatable :: wL(:),xS(:),dwL(:),dxS(:)
1175 double precision,
allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
1176 integer :: strtype,nstrb,nbb,nuni,nstr,bnx
1177 double precision :: qs,dxfirst,dxmid,lenstr
1179 integer :: iigrid,igrid,j,dir_loc
1180 double precision :: xbmin(1:ndim),xbmax(1:ndim)
1186 allocate(wl(numwl),dwl(numwl))
1189 wl(iwl)=wlmin+iwl*dwlg-half*dwlg
1221 call mpistop(
'Wrong direction_slit')
1224 allocate(xs(numxs),dxs(numxs),spectra(numwl,numxs),spectra_rc(numwl,numxs))
1226 allocate(wi(numwl,numxs,numwi))
1228 select case(strtype)
1230 dxs(:)=(xsmax-xsmin)/numxs
1232 xs(ixs)=xsmin+dxs(ixs)*(ixs-half)
1236 dxfirst=(xsmax-xsmin)*(one-qs)/(one-qs**numxs)
1239 dxs(ixs)=dxfirst*qs**(ixs-1)
1240 xs(ixs)=dxs(1)/(one-qs)*(one-qs**(ixs-1))+half*dxs(ixs)
1246 lenstr=(xsmax-xsmin)/(2.d0+nuni*(one-qs)/(one-qs**nstr))
1247 dxfirst=(xsmax-xsmin)/(dble(nuni)+2.d0/(one-qs)*(one-qs**nstr))
1253 dxfirst=lenstr*(one-qs)/(one-qs**nstr)
1256 if(nuni .gt. 0)
then
1257 do ixs=nstr+1,nstr+nuni
1259 xs(ixs)=lenstr+(dble(ixs)-0.5d0-nstr)*dxs(ixs)+xsmin
1264 dxs(ixs)=dxfirst*qs**(nstr-ixs)
1265 xs(ixs)=xsmin+lenstr-dxs(ixs)*half-dxfirst*(one-qs**(nstr-ixs))/(one-qs)
1268 do ixs=nstr+nuni+1,numxs
1269 dxs(ixs)=dxfirst*qs**(ixs-nstr-nuni-1)
1270 xs(ixs)=xsmax-lenstr+dxs(ixs)*half+dxfirst*(one-qs**(ixs-nstr-nuni-1))/(one-qs)
1273 call mpistop(
"unknown stretch type")
1296 call mpistop(
'Wrong combination of LOS and slit direction!')
1299 if (dir_loc==1)
then
1301 call mpistop(
'Wrong value for location_slit!')
1303 else if (dir_loc==2)
then
1305 call mpistop(
'Wrong value for location_slit!')
1309 call mpistop(
'Wrong value for location_slit!')
1315 do iigrid=1,igridstail; igrid=igrids(iigrid);
1325 call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
1329 if (spectra_rc(iwl,ixs)>smalldouble)
then
1330 wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
1337 call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
1339 deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
1346 integer,
intent(in) :: igrid,numWL,numXS,dir_loc
1348 double precision,
intent(in) :: wL(numWL),dwL(numWL)
1349 double precision,
intent(inout) :: spectra(numWL,numXS)
1351 integer :: direction_LOS
1352 integer :: ixO^L,ixI^L,ix^D,ixOnew
1353 double precision,
allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
1354 double precision :: wlc,wlwd
1357 double precision :: logTe,lineCent
1358 character (30) :: ion
1359 double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
1361 integer :: levelg,rft,ixSmin,ixSmax,iwL
1362 double precision :: flux_pix,dL
1374 ^d&ixomin^d=ixmlo^d\
1375 ^d&ixomax^d=ixmhi^d\
1376 ^d&iximin^d=
ixglo^d\
1377 ^d&iximax^d=
ixghi^d\
1378 allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
1381 if (dir_loc==1)
then
1382 do ix1=ixomin1,ixomax1
1390 else if (dir_loc==2)
then
1391 do ix2=ixomin2,ixomax2
1400 do ix3=ixomin3,ixomax3
1411 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
1412 v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
1413 call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
1414 if(
associated(fl%get_var_Rfactor))
then
1415 call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
1416 te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
1418 te(ixo^s)=pth(ixo^s)/(fl%Rfactor*rho(ixo^s))
1422 levelg=ps(igrid)%level
1425 {
do ix^d=ixomin^d,ixomax^d\}
1436 ixsmin=(block_nx1*(
node(pig1_,igrid)-1)+(ix1-ixomin1))*rft+1
1437 ixsmax=(block_nx1*(
node(pig1_,igrid)-1)+(ix1-ixomin1+1))*rft
1439 ixsmin=(block_nx2*(
node(pig2_,igrid)-1)+(ix2-ixomin2))*rft+1
1440 ixsmax=(block_nx2*(
node(pig2_,igrid)-1)+(ix2-ixomin2+1))*rft
1442 ixsmin=(block_nx3*(
node(pig3_,igrid)-1)+(ix3-ixomin3))*rft+1
1443 ixsmax=(block_nx3*(
node(pig3_,igrid)-1)+(ix3-ixomin3+1))*rft
1446 select case(direction_los)
1457 flux_pix=flux(ix^d)*wlrsl*dl*exp(-(wl(iwl)-wlc)**2/(2*wlwd**2))/(sqrt(2*
dpi)*wlwd)
1458 flux_pix=flux_pix*wslit/spacersl
1459 spectra(iwl,ixsmin:ixsmax)=spectra(iwl,ixsmin:ixsmax)+flux_pix
1464 deallocate(flux,v,pth,te,rho)
1471 integer,
intent(in) :: qunit
1473 character(20) :: datatype
1476 character (30) :: ion
1477 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1480 if (
mype==0) print *,
'###################################################'
1481 if (
mype==0) print *,
'Systhesizing EUV image'
1482 if (
mype==0)
write(*,
'(a,f8.3,a)')
' Wavelength: ',linecent,
' Angstrom'
1483 if (
mype==0) print *,
'Unit of EUV flux: DN s^-1 pixel^-1'
1484 if (
mype==0)
write(*,
'(a,f5.1,a,f5.1,a)')
' Pixel: ',spacersl*725.0,
' km x ',spacersl*725.0,
' km'
1486 '] of the simulation box is located at [X=0,Y=0] of the image'
1487 if (
mype==0) print *,
'Negative Doppler velocity indicates blueshift'
1489 datatype=
'image_euv'
1493 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
1495 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
1504 call mpistop(
'ERROR: Wrong LOS for synthesizing emission!')
1507 if (
mype==0) print *,
'Unit of length: arcsec (~725 km)'
1510 call mpistop(
'ERROR: Wrong resolution type')
1513 if (
mype==0) print *,
'###################################################'
1520 integer,
intent(in) :: qunit
1522 character(20) :: datatype
1524 if (
mype==0) print *,
'###################################################'
1525 if (
mype==0) print *,
'Systhesizing SXR image (observed at 1 AU).'
1527 if (
mype==0) print *,
'Unit of SXR flux: photons cm^-2 s^-1 pixel^-1'
1528 if (
mype==0)
write(*,
'(a,f6.1,a,f6.1,a)')
' Pixel: ',2.3*725.0,
' km x ',2.3*725.0,
' km'
1530 '] of the simulation box is located at [X=0,Y=0] of the image'
1532 datatype=
'image_sxr'
1536 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d3,
' km'
1538 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d5,
' km'
1547 call mpistop(
'ERROR: Wrong LOS for synthesizing emission!')
1550 if (
mype==0) print *,
'Unit of length: arcsec (~725 km)'
1553 call mpistop(
'ERROR: Wrong resolution type')
1556 if (
mype==0) print *,
'###################################################'
1566 integer,
intent(in) :: qunit
1568 character(20),
intent(in) :: datatype
1570 integer :: ix^D,numXI1,numXI2,numWI
1571 double precision :: xImin1,xImax1,xImin2,xImax2,xIcent1,xIcent2,dxI
1572 double precision,
allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:),wI(:,:,:)
1573 double precision,
allocatable :: EUVs(:,:),EUV(:,:),Dpls(:,:),Dpl(:,:)
1574 double precision,
allocatable :: SXRs(:,:),SXR(:,:)
1575 double precision :: vec_temp1(1:3),vec_temp2(1:3)
1576 double precision :: vec_z(1:3),vec_cor(1:3),xI_cor(1:2)
1577 double precision :: res,LOS_psi,r_max,r_loc
1580 character (30) :: ion
1581 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1582 double precision :: unitv,arcsec,RHESSI_rsl,pixel
1583 integer :: iigrid,igrid,i,j,numSI
1589 if (ix1==1) vec_cor(1)=xprobmin1
1590 if (ix1==2) vec_cor(1)=xprobmax1
1592 if (ix2==1) vec_cor(2)=xprobmin2
1593 if (ix2==2) vec_cor(2)=xprobmax2
1595 if (ix3==1) vec_cor(3)=xprobmin3
1596 if (ix3==2) vec_cor(3)=xprobmax3
1599 r_loc=r_loc+(vec_cor(2)-
x_origin(2))**2
1600 r_loc=r_loc+(vec_cor(3)-
x_origin(3))**2
1602 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1605 r_max=max(r_max,r_loc)
1609 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1615 ximin1=min(ximin1,xi_cor(1))
1616 ximax1=max(ximax1,xi_cor(1))
1617 ximin2=min(ximin2,xi_cor(2))
1618 ximax2=max(ximax2,xi_cor(2))
1630 xicent1=(ximin1+ximax1)/2.d0
1631 xicent2=(ximin2+ximax2)/2.d0
1639 if (datatype==
'image_euv')
then
1644 dxi=rhessi_rsl*arcsec
1646 numxi1=2*ceiling((ximax1-xicent1)/dxi/2.d0)
1647 ximin1=xicent1-numxi1*dxi
1648 ximax1=xicent1+numxi1*dxi
1650 numxi2=2*ceiling((ximax2-xicent2)/dxi/2.d0)
1651 ximin2=xicent2-numxi2*dxi
1652 ximax2=xicent2+numxi2*dxi
1654 allocate(xi1(numxi1),xi2(numxi2),dxi1(numxi1),dxi2(numxi2))
1656 xi1(ix1)=ximin1+dxi*(ix1-
half)
1660 xi2(ix2)=ximin2+dxi*(ix2-
half)
1665 if (datatype==
'image_euv')
then
1672 allocate(wi(numxi1,numxi2,numwi))
1673 allocate(euv(numxi1,numxi2),euvs(numxi1,numxi2))
1674 allocate(dpl(numxi1,numxi2),dpls(numxi1,numxi2))
1680 do iigrid=1,igridstail; igrid=igrids(iigrid);
1684 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision, &
1686 call mpi_allreduce(dpls,dpl,numsi,mpi_double_precision, &
1691 if(euv(ix1,ix2)/=0)
then
1692 dpl(ix1,ix2)=(dpl(ix1,ix2)/euv(ix1,ix2))*unitv
1697 wi(ix1,ix2,1)=euv(ix1,ix2)
1698 wi(ix1,ix2,2)=dpl(ix1,ix2)
1706 call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
1707 deallocate(wi,euv,euvs,dpl,dpls)
1710 if (datatype==
'image_sxr')
then
1712 allocate(wi(numxi1,numxi2,numwi))
1713 allocate(sxr(numxi1,numxi2),sxrs(numxi1,numxi2))
1717 do iigrid=1,igridstail; igrid=igrids(iigrid);
1721 call mpi_allreduce(sxrs,sxr,numsi,mpi_double_precision, &
1735 call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
1736 deallocate(wi,sxr,sxrs)
1739 deallocate(xi1,xi2,dxi1,dxi2)
1744 integer,
intent(in) :: igrid,numXI1,numXI2
1745 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
1746 double precision,
intent(in) :: dxI
1748 double precision,
intent(out) :: SXR(numXI1,numXI2)
1750 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
1751 double precision :: xb^L,xd^D
1752 double precision,
allocatable :: flux(:^D&)
1753 double precision :: vloc(1:3),res
1754 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
1755 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC
1756 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2)
1758 double precision :: sigma_PSF,RHESSI_rsl,sigma0,factor
1759 double precision :: arcsec,pixel,area_1AU
1761 ^d&ixomin^d=ixmlo^d\
1762 ^d&ixomax^d=ixmhi^d\
1763 ^d&iximin^d=
ixglo^d\
1764 ^d&iximax^d=
ixghi^d\
1768 allocate(flux(ixi^s))
1780 pixel=rhessi_rsl*arcsec
1781 sigma0=sigma_psf*pixel
1783 {
do ix^d=ixomin^d,ixomax^d\}
1785 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi1(^d))/(dxi/4.d0)));
1786 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi2(^d))/(dxi/4.d0)));
1787 ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
1789 {
do isubc^d=1,nsubc^d\}
1790 ^d&xsubc(^d)=ps(igrid)%x(ix^dd,^d)-half*ps(igrid)%dx(ix^dd,^d)+(isubc^d-half)*dxsubc^d;
1791 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*
unit_length**3/area_1au
1794 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
1795 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
1796 ixpmin1=max(1,ixp1-3)
1797 ixpmax1=min(ixp1+3,numxi1)
1798 ixpmin2=max(1,ixp2-3)
1799 ixpmax2=min(ixp2+3,numxi2)
1800 do ixp1=ixpmin1,ixpmax1
1801 do ixp2=ixpmin2,ixpmax2
1802 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
1803 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
1804 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
1805 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
1806 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
1807 sxr(ixp1,ixp2)=sxr(ixp1,ixp2)+fluxsubc*factor
1817 integer,
intent(in) :: igrid,numXI1,numXI2
1818 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
1819 double precision,
intent(in) :: dxI
1821 double precision,
intent(out) :: EUV(numXI1,numXI2),Dpl(numXI1,numXI2)
1823 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
1824 double precision :: xb^L,xd^D
1825 double precision,
allocatable :: flux(:^D&),v(:^D&),rho(:^D&)
1826 double precision :: vloc(1:3),res
1827 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
1828 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC
1829 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2)
1832 double precision :: logTe
1833 character (30) :: ion
1834 double precision :: lineCent
1835 double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
1836 double precision :: unitv,arcsec,pixel
1838 ^d&ixomin^d=ixmlo^d\
1839 ^d&ixomax^d=ixmhi^d\
1840 ^d&iximin^d=
ixglo^d\
1841 ^d&iximax^d=
ixghi^d\
1845 allocate(flux(ixi^s),v(ixi^s),rho(ixi^s))
1848 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
1849 {
do ix^d=ixomin^d,ixomax^d\}
1851 vloc(j)=ps(igrid)%w(ix^d,iw_mom(j))/rho(ix^d)
1865 pixel=spacersl*arcsec
1866 sigma0=sigma_psf*pixel
1867 {
do ix^d=ixomin^d,ixomax^d\}
1869 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi1(^d))/(dxi/2.d0)));
1870 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi2(^d))/(dxi/2.d0)));
1871 ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
1873 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*
unit_length*1.d2/dxi/dxi
1875 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*
unit_length/dxi/dxi
1878 {
do isubc^d=1,nsubc^d\}
1879 ^d&xsubc(^d)=ps(igrid)%x(ix^dd,^d)-half*ps(igrid)%dx(ix^dd,^d)+(isubc^d-half)*dxsubc^d;
1883 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
1884 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
1885 ixpmin1=max(1,ixp1-3)
1886 ixpmax1=min(ixp1+3,numxi1)
1887 ixpmin2=max(1,ixp2-3)
1888 ixpmax2=min(ixp2+3,numxi2)
1889 do ixp1=ixpmin1,ixpmax1
1890 do ixp2=ixpmin2,ixpmax2
1891 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
1892 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
1893 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
1894 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
1895 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
1896 euv(ixp1,ixp2)=euv(ixp1,ixp2)+fluxsubc*factor
1897 dpl(ixp1,ixp2)=dpl(ixp1,ixp2)+fluxsubc*factor*v(ix^d)
1912 integer,
intent(in) :: qunit
1913 character(20),
intent(in) :: datatype
1916 double precision :: dx^D
1917 integer :: numX^D,ix^D
1918 double precision,
allocatable :: EUV(:,:),EUVs(:,:),Dpl(:,:),Dpls(:,:)
1919 double precision,
allocatable :: SXR(:,:),SXRs(:,:),wI(:,:,:)
1920 double precision,
allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:),dxIi
1921 integer :: numXI1,numXI2,numSI,numWI
1922 double precision :: xI^L
1923 integer :: iigrid,igrid,i,j
1924 double precision,
allocatable :: xIF1(:),xIF2(:),dxIF1(:),dxIF2(:)
1925 integer :: nXIF1,nXIF2
1926 double precision :: xIF^L
1928 double precision :: unitv,arcsec,RHESSI_rsl
1929 integer :: strtype^D,nstrb^D,nbb^D,nuni^D,nstr^D,bnx^D
1930 double precision :: qs^D,dxfirst^D,dxmid^D,lenstr^D
1954 if (
mype==0)
write(*,
'(a)')
' LOS vector: [-1.00 0.00 0.00]'
1955 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [ 0.00 1.00 0.00]'
1956 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 0.00 1.00]'
1974 if (
mype==0)
write(*,
'(a)')
' LOS vector: [ 0.00 -1.00 0.00]'
1975 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [-1.00 0.00 0.00]'
1976 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 0.00 1.00]'
1994 if (
mype==0)
write(*,
'(a)')
' LOS vector: [ 0.00 0.00 -1.00]'
1995 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [ 1.00 0.00 0.00]'
1996 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 1.00 0.00]'
1998 allocate(xif1(nxif1),xif2(nxif2),dxif1(nxif1),dxif2(nxif2))
2001 select case(strtype1)
2003 dxif1(:)=(xifmax1-xifmin1)/nxif1
2005 xif1(ix1)=xifmin1+dxif1(ix1)*(ix1-
half)
2009 dxfirst1=(xifmax1-xifmin1)*(
one-qs1)/(
one-qs1**nxif1)
2012 dxif1(ix1)=dxfirst1*qs1**(ix1-1)
2013 xif1(ix1)=dxif1(1)/(
one-qs1)*(
one-qs1**(ix1-1))+
half*dxif1(ix1)
2018 nuni1=nbb1-nstrb1*bnx1
2019 lenstr1=(xifmax1-xifmin1)/(2.d0+nuni1*(
one-qs1)/(
one-qs1**nstr1))
2020 dxfirst1=(xifmax1-xifmin1)/(dble(nuni1)+2.d0/(
one-qs1)*(
one-qs1**nstr1))
2026 dxfirst1=lenstr1*(
one-qs1)/(
one-qs1**nstr1)
2029 if(nuni1 .gt. 0)
then
2030 do ix1=nstr1+1,nstr1+nuni1
2032 xif1(ix1)=lenstr1+(dble(ix1)-0.5d0-nstr1)*dxif1(ix1)+xifmin1
2037 dxif1(ix1)=dxfirst1*qs1**(nstr1-ix1)
2038 xif1(ix1)=xifmin1+lenstr1-dxif1(ix1)*
half-dxfirst1*(
one-qs1**(nstr1-ix1))/(
one-qs1)
2041 do ix1=nstr1+nuni1+1,nxif1
2042 dxif1(ix1)=dxfirst1*qs1**(ix1-nstr1-nuni1-1)
2043 xif1(ix1)=xifmax1-lenstr1+dxif1(ix1)*
half+dxfirst1*(
one-qs1**(ix1-nstr1-nuni1-1))/(
one-qs1)
2046 call mpistop(
"unknown stretch type")
2049 select case(strtype2)
2051 dxif2(:)=(xifmax2-xifmin2)/nxif2
2053 xif2(ix2)=xifmin2+dxif2(ix2)*(ix2-
half)
2057 dxfirst2=(xifmax2-xifmin2)*(
one-qs2)/(
one-qs2**nxif2)
2060 dxif2(ix2)=dxfirst2*qs2**(ix2-1)
2061 xif2(ix2)=dxif2(1)/(
one-qs2)*(
one-qs2**(ix2-1))+
half*dxif2(ix2)
2066 nuni2=nbb2-nstrb2*bnx2
2067 lenstr2=(xifmax2-xifmin2)/(2.d0+nuni2*(
one-qs2)/(
one-qs2**nstr2))
2068 dxfirst2=(xifmax2-xifmin2)/(dble(nuni2)+2.d0/(
one-qs2)*(
one-qs2**nstr2))
2074 dxfirst2=lenstr2*(
one-qs2)/(
one-qs2**nstr2)
2077 if(nuni2 .gt. 0)
then
2078 do ix2=nstr2+1,nstr2+nuni2
2080 xif2(ix2)=lenstr2+(dble(ix2)-0.5d0-nstr2)*dxif2(ix2)+xifmin2
2085 dxif2(ix2)=dxfirst2*qs2**(nstr2-ix2)
2086 xif2(ix2)=xifmin2+lenstr2-dxif2(ix2)*
half-dxfirst2*(
one-qs2**(nstr2-ix2))/(
one-qs2)
2089 do ix2=nstr2+nuni2+1,nxif2
2090 dxif2(ix2)=dxfirst2*qs2**(ix2-nstr2-nuni2-1)
2091 xif2(ix2)=xifmax2-lenstr2+dxif2(ix2)*
half+dxfirst2*(
one-qs2**(ix2-nstr2-nuni2-1))/(
one-qs2)
2094 call mpistop(
"unknown stretch type")
2098 if (datatype==
'image_euv')
then
2105 allocate(wi(nxif1,nxif2,numwi))
2106 allocate(euvs(nxif1,nxif2),euv(nxif1,nxif2))
2107 allocate(dpl(nxif1,nxif2),dpls(nxif1,nxif2))
2112 do iigrid=1,igridstail; igrid=igrids(iigrid);
2116 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision, &
2118 call mpi_allreduce(dpls,dpl,numsi,mpi_double_precision, &
2123 if(euv(ix1,ix2)/=0)
then
2124 dpl(ix1,ix2)=(dpl(ix1,ix2)/euv(ix1,ix2))*unitv
2134 call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
2135 deallocate(wi,euv,euvs,dpl,dpls)
2139 if (datatype==
'image_sxr')
then
2147 allocate(wi(nxif1,nxif2,numwi))
2148 allocate(sxrs(nxif1,nxif2),sxr(nxif1,nxif2))
2151 do iigrid=1,igridstail; igrid=igrids(iigrid);
2155 call mpi_allreduce(sxrs,sxr,numsi,mpi_double_precision, &
2158 sxr=sxr*(rhessi_rsl*arcsec)**2
2166 call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
2167 deallocate(wi,sxr,sxrs)
2170 deallocate(xif1,xif2,dxif1,dxif2)
2177 integer,
intent(in) :: igrid,nXIF1,nXIF2
2178 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
2179 double precision,
intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
2181 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
2183 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
2184 double precision :: xb^L,xd^D
2185 double precision,
allocatable :: flux(:^D&),v(:^D&),rho(:^D&)
2186 double precision,
allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
2187 double precision,
allocatable :: EUVg(:,:),Fvg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
2188 integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
2189 double precision :: EUVt,Fvt,xc^L,xg^L,r2
2190 integer :: ixP^L,ixP^D
2191 integer :: direction_LOS
2201 ^d&ixomin^d=ixmlo^d\
2202 ^d&ixomax^d=ixmhi^d\
2203 ^d&iximin^d=
ixglo^d\
2204 ^d&iximax^d=
ixghi^d\
2208 allocate(flux(ixi^s),v(ixi^s),rho(ixi^s))
2209 allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
2210 dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
2211 dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
2212 dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
2215 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
2216 v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
2220 levelg=ps(igrid)%level
2224 select case(direction_los)
2235 allocate(euvg(nxg1,nxg2),fvg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
2242 select case(direction_los)
2244 do ix2=ixomin2,ixomax2
2245 ixgmin1=(ix2-1)*rft+1
2247 do ix3=ixomin3,ixomax3
2248 ixgmin2=(ix3-1)*rft+1
2252 do ix1=ixomin1,ixomax1
2256 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
2257 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
2261 do ix3=ixomin3,ixomax3
2262 ixgmin1=(ix3-1)*rft+1
2264 do ix1=ixomin1,ixomax1
2265 ixgmin2=(ix1-1)*rft+1
2269 do ix2=ixomin2,ixomax2
2273 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
2274 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
2278 do ix1=ixomin1,ixomax1
2279 ixgmin1=(ix1-1)*rft+1
2281 do ix2=ixomin2,ixomax2
2282 ixgmin2=(ix2-1)*rft+1
2286 do ix3=ixomin3,ixomax3
2290 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
2291 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
2302 select case(direction_los)
2304 ixgmin1=(ixomin2-1)*rft+1
2306 ixgmin2=(ixomin3-1)*rft+1
2309 ixgmin1=(ixomin3-1)*rft+1
2311 ixgmin2=(ixomin1-1)*rft+1
2314 ixgmin1=(ixomin1-1)*rft+1
2316 ixgmin2=(ixomin2-1)*rft+1
2320 select case(direction_los)
2322 ixpmin1=(
node(pig2_,igrid)-1)*rft*block_nx2+1
2323 ixpmax1=
node(pig2_,igrid)*rft*block_nx2
2324 ixpmin2=(
node(pig3_,igrid)-1)*rft*block_nx3+1
2325 ixpmax2=
node(pig3_,igrid)*rft*block_nx3
2327 ixpmin1=(
node(pig3_,igrid)-1)*rft*block_nx3+1
2328 ixpmax1=
node(pig3_,igrid)*rft*block_nx3
2329 ixpmin2=(
node(pig1_,igrid)-1)*rft*block_nx1+1
2330 ixpmax2=
node(pig1_,igrid)*rft*block_nx1
2332 ixpmin1=(
node(pig1_,igrid)-1)*rft*block_nx1+1
2333 ixpmax1=
node(pig1_,igrid)*rft*block_nx1
2334 ixpmin2=(
node(pig2_,igrid)-1)*rft*block_nx2+1
2335 ixpmax2=
node(pig2_,igrid)*rft*block_nx2
2337 xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
2338 xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
2339 dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
2340 dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
2341 euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
2342 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
2343 dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
2344 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
2346 deallocate(flux,v,dxb1,dxb2,dxb3,euvg,fvg,xg1,xg2,dxg1,dxg2)
2353 integer,
intent(in) :: igrid,nXIF1,nXIF2
2354 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
2355 double precision,
intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
2357 double precision,
intent(out) :: SXR(nXIF1,nXIF2)
2359 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
2360 double precision :: xb^L,xd^D
2361 double precision,
allocatable :: flux(:^D&)
2362 double precision,
allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
2363 double precision,
allocatable :: SXRg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
2364 integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
2365 double precision :: SXRt,xc^L,xg^L,r2,area_1AU
2366 integer :: ixP^L,ixP^D
2367 integer :: direction_LOS
2377 ^d&ixomin^d=ixmlo^d\
2378 ^d&ixomax^d=ixmhi^d\
2379 ^d&iximin^d=
ixglo^d\
2380 ^d&iximax^d=
ixghi^d\
2384 allocate(flux(ixi^s))
2385 allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
2386 dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
2387 dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
2388 dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
2393 levelg=ps(igrid)%level
2397 select case(direction_los)
2408 allocate(sxrg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
2414 select case(direction_los)
2416 do ix2=ixomin2,ixomax2
2417 ixgmin1=(ix2-1)*rft+1
2419 do ix3=ixomin3,ixomax3
2420 ixgmin2=(ix3-1)*rft+1
2423 do ix1=ixomin1,ixomax1
2426 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
2430 do ix3=ixomin3,ixomax3
2431 ixgmin1=(ix3-1)*rft+1
2433 do ix1=ixomin1,ixomax1
2434 ixgmin2=(ix1-1)*rft+1
2437 do ix2=ixomin2,ixomax2
2440 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
2444 do ix1=ixomin1,ixomax1
2445 ixgmin1=(ix1-1)*rft+1
2447 do ix2=ixomin2,ixomax2
2448 ixgmin2=(ix2-1)*rft+1
2451 do ix3=ixomin3,ixomax3
2454 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
2464 select case(direction_los)
2466 ixgmin1=(ixomin2-1)*rft+1
2468 ixgmin2=(ixomin3-1)*rft+1
2471 ixgmin1=(ixomin3-1)*rft+1
2473 ixgmin2=(ixomin1-1)*rft+1
2476 ixgmin1=(ixomin1-1)*rft+1
2478 ixgmin2=(ixomin2-1)*rft+1
2482 select case(direction_los)
2484 ixpmin1=(
node(pig2_,igrid)-1)*rft*block_nx2+1
2485 ixpmax1=
node(pig2_,igrid)*rft*block_nx2
2486 ixpmin2=(
node(pig3_,igrid)-1)*rft*block_nx3+1
2487 ixpmax2=
node(pig3_,igrid)*rft*block_nx3
2489 ixpmin1=(
node(pig3_,igrid)-1)*rft*block_nx3+1
2490 ixpmax1=
node(pig3_,igrid)*rft*block_nx3
2491 ixpmin2=(
node(pig1_,igrid)-1)*rft*block_nx1+1
2492 ixpmax2=
node(pig1_,igrid)*rft*block_nx1
2494 ixpmin1=(
node(pig1_,igrid)-1)*rft*block_nx1+1
2495 ixpmax1=
node(pig1_,igrid)*rft*block_nx1
2496 ixpmin2=(
node(pig2_,igrid)-1)*rft*block_nx2+1
2497 ixpmax2=
node(pig2_,igrid)*rft*block_nx2
2499 xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
2500 xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
2501 dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
2502 dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
2503 sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
2504 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
2506 deallocate(flux,dxb1,dxb2,dxb3,sxrg,xg1,xg2,dxg1,dxg2)
2510 subroutine output_data(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,datatype)
2514 integer,
intent(in) :: qunit,nXO1,nXO2,nWO
2515 double precision,
intent(in) :: dxO1(nxO1),dxO2(nxO2)
2516 double precision,
intent(in) :: xO1(nXO1),xO2(nxO2)
2517 double precision,
intent(inout) :: wO(nXO1,nXO2,nWO)
2518 character(20),
intent(in) :: datatype
2520 integer :: nPiece,nP1,nP2,nC1,nC2,nWC
2521 integer :: piece_nmax1,piece_nmax2,ix1,ix2,j,ipc,ixc1,ixc2
2522 double precision,
allocatable :: xC(:,:,:,:),wC(:,:,:,:),dxC(:,:,:,:)
2523 character(len=std_len) :: resolution_type
2529 if (abs(wo(ix1,ix2,j))<smalldouble) wo(ix1,ix2,j)=zero
2537 piece_nmax1=block_nx2
2538 piece_nmax2=block_nx3
2540 piece_nmax1=block_nx3
2541 piece_nmax2=block_nx1
2543 piece_nmax1=block_nx1
2544 piece_nmax2=block_nx2
2546 else if (datatype==
'image_sxr' .and.
resolution_sxr==
'data')
then
2548 piece_nmax1=block_nx2
2549 piece_nmax2=block_nx3
2551 piece_nmax1=block_nx3
2552 piece_nmax2=block_nx1
2554 piece_nmax1=block_nx1
2555 piece_nmax2=block_nx2
2560 piece_nmax2=block_nx1
2562 piece_nmax2=block_nx2
2564 piece_nmax2=block_nx3
2571 loopn1:
do j=piece_nmax1,1,-1
2572 if(mod(nxo1,j)==0)
then
2577 loopn2:
do j=piece_nmax2,1,-1
2578 if(mod(nxo2,j)==0)
then
2590 case(
'EIvtuCCmpi',
'ESvtuCCmpi',
'SIvtuCCmpi')
2592 allocate(xc(npiece,nc1,nc2,2))
2593 allocate(dxc(npiece,nc1,nc2,2))
2594 allocate(wc(npiece,nc1,nc2,nwo))
2598 ix1=mod(ipc-1,np1)*nc1+ixc1
2599 ix2=floor(1.0*(ipc-1)/np1)*nc2+ixc2
2600 xc(ipc,ixc1,ixc2,1)=xo1(ix1)
2601 xc(ipc,ixc1,ixc2,2)=xo2(ix2)
2602 dxc(ipc,ixc1,ixc2,1)=dxo1(ix1)
2603 dxc(ipc,ixc1,ixc2,2)=dxo2(ix2)
2605 wc(ipc,ixc1,ixc2,j)=wo(ix1,ix2,j)
2612 deallocate(xc,dxc,wc)
2613 case(
'EIvtiCCmpi',
'ESvtiCCmpi',
'SIvtiCCmpi')
2617 if (sum(
stretch_type(:))>0 .and. resolution_type==
'data')
then
2618 call mpistop(
"Error in synthesize emission: vti is not supported for data resolution")
2620 call write_image_vticc(qunit,xo1,xo2,dxo1,dxo2,wo,nxo1,nxo2,nwo,nc1,nc2)
2623 call mpistop(
"Error in synthesize emission: Unknown convert_type")
2629 subroutine write_image_vticc(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,nC1,nC2)
2633 integer,
intent(in) :: qunit,nXO1,nXO2,nWO,nC1,nC2
2634 double precision,
intent(in) :: xO1(nXO1),xO2(nxO2)
2635 double precision,
intent(in) :: dxO1(nxO1),dxO2(nxO2)
2636 double precision,
intent(in) :: wO(nXO1,nXO2,nWO)
2638 double precision :: origin(1:3), spacing(1:3)
2639 integer :: wholeExtent(1:6),extent(1:6)
2640 integer :: nP1,nP2,iP1,iP2,iw
2641 integer :: ixC1,ixC2,ixCmin1,ixCmax1,ixCmin2,ixCmax2
2645 character (70) :: subname,wname,vname,nameL,nameS
2646 character (len=std_len) :: filename
2648 double precision :: logTe
2649 character (30) :: ion
2650 double precision :: line_center
2651 double precision :: spatial_rsl,spectral_rsl,sigma_PSF,wslit
2654 origin(1)=xo1(1)-0.5d0*dxo1(1)
2655 origin(2)=xo2(1)-0.5d0*dxo2(1)
2674 inquire(qunit,opened=fileopen)
2675 if(.not.fileopen)
then
2680 write(filename,
'(a,i4.4,a)') trim(
filename_euv),filenr,
".vti"
2682 write(filename,
'(a,i4.4,a)') trim(
filename_sxr),filenr,
".vti"
2686 open(qunit,file=filename,status=
'unknown',form=
'formatted')
2690 write(qunit,
'(a)')
'<?xml version="1.0"?>'
2691 write(qunit,
'(a)',advance=
'no')
'<VTKFile type="ImageData"'
2692 write(qunit,
'(a)')
' version="0.1" byte_order="LittleEndian">'
2693 write(qunit,
'(a,3(1pe14.6),a,6(i10),a,3(1pe14.6),a)')
' <ImageData Origin="',&
2694 origin,
'" WholeExtent="',wholeextent,
'" Spacing="',spacing,
'">'
2696 write(qunit,
'(a)')
'<FieldData>'
2697 write(qunit,
'(2a)')
'<DataArray type="Float32" Name="TIME" ',&
2698 'NumberOfTuples="1" format="ascii">'
2700 write(qunit,
'(a)')
'</DataArray>'
2702 write(qunit,
'(2a)')
'<DataArray type="Float32" Name="logT" ',&
2703 'NumberOfTuples="1" format="ascii">'
2704 write(qunit,*) real(logte)
2705 write(qunit,
'(a)')
'</DataArray>'
2707 write(qunit,
'(a)')
'</FieldData>'
2712 extent(1)=(ip1-1)*nc1
2714 extent(3)=(ip2-1)*nc2
2720 write(qunit,
'(a,6(i10),a)') &
2721 '<Piece Extent="',extent,
'">'
2722 write(qunit,
'(a)')
'<CellData>'
2735 if (iw==2) vname=
'Doppler_velocity'
2751 write(qunit,
'(a,a,a)')&
2752 '<DataArray type="Float64" Name="',trim(vname),
'" format="ascii">'
2753 write(qunit,
'(200(1pe14.6))') ((wo(ixc1,ixc2,iw),ixc1=ixcmin1,ixcmax1),ixc2=ixcmin2,ixcmax2)
2754 write(qunit,
'(a)')
'</DataArray>'
2756 write(qunit,
'(a)')
'</CellData>'
2757 write(qunit,
'(a)')
'</Piece>'
2761 write(qunit,
'(a)')
'</ImageData>'
2762 write(qunit,
'(a)')
'</VTKFile>'
2772 integer,
intent(in) :: qunit
2773 integer,
intent(in) :: nPiece,nC1,nC2,nWC
2774 double precision,
intent(in) :: xC(nPiece,nC1,nC2,2),dxC(nPiece,nc1,nc2,2)
2775 double precision,
intent(in) :: wC(nPiece,nC1,nC2,nWC)
2776 character(20),
intent(in) :: datatype
2779 double precision :: xP(nPiece,nC1+1,nC2+1,2)
2782 character (70) :: subname,wname,vname,nameL,nameS
2783 character (len=std_len) :: filename
2784 integer :: ixC1,ixC2,ixP,ix1,ix2,j
2785 integer :: nc,np,icel,VTK_type
2787 double precision :: logTe
2788 character (30) :: ion
2789 double precision :: line_center
2790 double precision :: spatial_rsl,spectral_rsl,sigma_PSF,wslit
2800 if (ix1<np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1,1,1)-0.5d0*dxc(ixp,ix1,1,1)
2801 if (ix1==np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1-1,1,1)+0.5d0*dxc(ixp,ix1-1,1,1)
2802 if (ix2<np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2,2)-0.5d0*dxc(ixp,1,ix2,2)
2803 if (ix2==np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2-1,2)+0.5d0*dxc(ixp,1,ix2-1,2)
2808 if (datatype==
'image_euv')
then
2810 else if (datatype==
'spectrum_euv')
then
2815 inquire(qunit,opened=fileopen)
2816 if(.not.fileopen)
then
2820 if (datatype==
'image_euv')
then
2821 write(filename,
'(a,i4.4,a)') trim(
filename_euv),filenr,
".vtu"
2822 else if (datatype==
'image_sxr')
then
2823 write(filename,
'(a,i4.4,a)') trim(
filename_sxr),filenr,
".vtu"
2824 else if (datatype==
'spectrum_euv')
then
2827 open(qunit,file=filename,status=
'unknown',form=
'formatted')
2830 write(qunit,
'(a)')
'<?xml version="1.0"?>'
2831 write(qunit,
'(a)',advance=
'no')
'<VTKFile type="UnstructuredGrid"'
2832 write(qunit,
'(a)')
' version="0.1" byte_order="LittleEndian">'
2833 write(qunit,
'(a)')
'<UnstructuredGrid>'
2834 write(qunit,
'(a)')
'<FieldData>'
2835 write(qunit,
'(2a)')
'<DataArray type="Float32" Name="TIME" ',&
2836 'NumberOfTuples="1" format="ascii">'
2838 write(qunit,
'(a)')
'</DataArray>'
2839 if (datatype==
'image_euv' .or. datatype==
'spectrum_euv')
then
2840 write(qunit,
'(2a)')
'<DataArray type="Float32" Name="logT" ',&
2841 'NumberOfTuples="1" format="ascii">'
2842 write(qunit,*) real(logte)
2843 write(qunit,
'(a)')
'</DataArray>'
2845 write(qunit,
'(a)')
'</FieldData>'
2847 write(qunit,
'(a,i7,a,i7,a)') &
2848 '<Piece NumberOfPoints="',np,
'" NumberOfCells="',nc,
'">'
2849 write(qunit,
'(a)')
'<CellData>'
2851 if (datatype==
'image_euv')
then
2861 if (j==2) vname=
'Doppler_velocity'
2862 else if (datatype==
'image_sxr')
then
2870 else if (datatype==
'spectrum_euv')
then
2877 write(qunit,
'(a,a,a)')&
2878 '<DataArray type="Float64" Name="',trim(vname),
'" format="ascii">'
2879 write(qunit,
'(200(1pe14.6))') ((wc(ixp,ixc1,ixc2,j),ixc1=1,nc1),ixc2=1,nc2)
2880 write(qunit,
'(a)')
'</DataArray>'
2882 write(qunit,
'(a)')
'</CellData>'
2883 write(qunit,
'(a)')
'<Points>'
2884 write(qunit,
'(a)')
'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
2889 write(qunit,
'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
2891 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
2893 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
2895 else if (datatype==
'image_sxr' .and.
resolution_sxr==
'data')
then
2897 write(qunit,
'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
2899 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
2901 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
2904 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
2908 write(qunit,
'(a)')
'</DataArray>'
2909 write(qunit,
'(a)')
'</Points>'
2911 write(qunit,
'(a)')
'<Cells>'
2912 write(qunit,
'(a)')
'<DataArray type="Int32" Name="connectivity" format="ascii">'
2915 write(qunit,
'(4(i7))') ix1-1+(ix2-1)*np1,ix1+(ix2-1)*np1,&
2916 ix1-1+ix2*np1,ix1+ix2*np1
2919 write(qunit,
'(a)')
'</DataArray>'
2921 write(qunit,
'(a)')
'<DataArray type="Int32" Name="offsets" format="ascii">'
2923 write(qunit,
'(i7)') icel*(2**2)
2925 write(qunit,
'(a)')
'</DataArray>'
2927 write(qunit,
'(a)')
'<DataArray type="Int32" Name="types" format="ascii">'
2931 write(qunit,
'(i2)') vtk_type
2933 write(qunit,
'(a)')
'</DataArray>'
2934 write(qunit,
'(a)')
'</Cells>'
2935 write(qunit,
'(a)')
'</Piece>'
2937 write(qunit,
'(a)')
'</UnstructuredGrid>'
2938 write(qunit,
'(a)')
'</VTKFile>'
2944 double precision,
intent(in) :: vec_in1(1:3),vec_in2(1:3)
2945 double precision,
intent(out) :: res
2950 res=res+vec_in1(j)*vec_in2(j)
2956 double precision,
intent(in) :: vec_in1(1:3),vec_in2(1:3)
2957 double precision,
intent(out) :: vec_out(1:3)
2959 vec_out(1)=vec_in1(2)*vec_in2(3)-vec_in1(3)*vec_in2(2)
2960 vec_out(2)=vec_in1(3)*vec_in2(1)-vec_in1(1)*vec_in2(3)
2961 vec_out(3)=vec_in1(1)*vec_in2(2)-vec_in1(2)*vec_in2(1)
2967 double precision :: LOS_psi
2968 double precision :: vec_z(1:3),vec_temp1(1:3),vec_temp2(1:3)
2991 vec_xi1=vec_temp1*cos(los_psi)-vec_temp2*sin(los_psi)
2992 vec_xi2=vec_temp2*cos(los_psi)+vec_temp1*sin(los_psi)
3006 double precision :: x_3D(1:3),x_image(1:2)
3007 double precision :: res,res_origin
3011 x_image(1)=res-res_origin
3014 x_image(2)=res-res_origin
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
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
This module contains definitions of global parameters and variables and some generic functions/subrou...
character(len=std_len) filename_sxr
Base file name for synthetic SXR emission output.
integer spectrum_wl
wave length for spectrum
integer ixghi
Upper index of grid block arrays.
character(len=std_len) filename_spectrum
Base file name for synthetic EUV spectrum output.
double precision global_time
The global simulation time.
integer snapshotini
Resume from the snapshot with this index.
character(len=std_len) filename_euv
Base file name for synthetic EUV emission output.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) 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.
integer mype
The rank of the current MPI task.
integer, dimension(:), allocatable, parameter d
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
double precision image_rotate
rotation of image
character(len=std_len) resolution_sxr
resolution of the SXR image
integer snapshotnext
IO: snapshot and collapsed views output numbers/labels.
integer, dimension(ndim) nstretchedblocks_baselevel
(even) number of (symmetrically) stretched blocks at AMR level 1, per dimension
double precision unit_velocity
Physical scaling factor for velocity.
double precision, dimension(ndim) qstretch_baselevel
stretch factor between cells at AMR level 1, per dimension
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision los_theta
direction of the line of sight (LOS)
double precision spectrum_window_max
character(len=std_len) resolution_euv
resolution of the EUV image
integer wavelength
wavelength for output
integer, dimension(ndim) stretch_type
What kind of stretching is used per dimension.
integer, parameter rpxmax
character(len=std_len) resolution_spectrum
resolution of the spectrum
logical big_image
big image
double precision spectrum_window_min
spectral window
integer refine_max_level
Maximal number of AMR levels.
integer direction_slit
direction of the slit
double precision, dimension(1:3) x_origin
where the is the origin (X=0,Y=0) of image
integer, dimension(:,:), allocatable node
double precision, dimension(ndim) dxlevel
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 integrate_euv_data_resol(igrid, nXIF1, nXIF2, xIF1, xIF2, dxIF1, dxIF2, fl, EUV, Dpl)
double precision, dimension(1:101) f_304
double precision, dimension(1:101) f_193
double precision, dimension(1:60) f_264
double precision, dimension(1:60) f_263
subroutine get_euv_image(qunit, fl)
subroutine get_line_info(wl, ion, mass, logTe, line_center, spatial_px, spectral_px, sigma_PSF, width_slit)
double precision, dimension(1:60) t_eis1
subroutine integrate_spectra_data_resol(igrid, wL, dwL, spectra, numWL, numXS, dir_loc, fl)
double precision, dimension(1:101) f_131
subroutine get_image_inst_resol(qunit, datatype, fl)
double precision, dimension(1:60) f_255
subroutine integrate_sxr_inst_resol(igrid, numXI1, numXI2, xI1, xI2, dxI, fl, SXR)
subroutine write_image_vtucc(qunit, xC, wC, dxC, nPiece, nC1, nC2, nWC, datatype)
subroutine get_euv(wl, ixIL, ixOL, w, x, fl, flux)
subroutine get_cor_image(x_3D, x_image)
subroutine write_image_vticc(qunit, xO1, xO2, dxO1, dxO2, wO, nXO1, nXO2, nWO, nC1, nC2)
double precision, dimension(1:101) t_aia
subroutine get_spectrum_inst_resol(qunit, datatype, fl)
subroutine integrate_euv_inst_resol(igrid, numXI1, numXI2, xI1, xI2, dxI, fl, EUV, Dpl)
double precision, dimension(1:60) t_eis2
subroutine get_spectrum_data_resol(qunit, datatype, fl)
double precision, dimension(1:101) f_171
subroutine output_data(qunit, xO1, xO2, dxO1, dxO2, wO, nXO1, nXO2, nWO, datatype)
double precision, dimension(1:101) f_94
double precision, dimension(1:3) vec_xi1
subroutine get_image_data_resol(qunit, datatype, fl)
subroutine dot_product_loc(vec_in1, vec_in2, res)
subroutine get_sxr_image(qunit, fl)
double precision, dimension(1:60) f_192
subroutine integrate_sxr_data_resol(igrid, nXIF1, nXIF2, xIF1, xIF2, dxIF1, dxIF2, fl, SXR)
double precision, dimension(1:3) vec_xi2
subroutine init_vectors()
double precision, dimension(1:41) f_1354
subroutine cross_product_loc(vec_in1, vec_in2, vec_out)
double precision, dimension(1:101) f_335
subroutine get_sxr(ixIL, ixOL, w, x, fl, flux, El, Eu)
double precision, dimension(1:41) t_iris
subroutine get_euv_spectrum(qunit, fl)
double precision, dimension(1:101) f_211
subroutine integrate_spectra_inst_resol(igrid, wL, dwLg, xS, dxSg, spectra, numWL, numXS, fl)
subroutine get_goes_flux_grid(ixIL, ixOL, w, x, dV, xboxL, xbL, fl, eflux_grid)
subroutine get_goes_sxr_flux(xboxL, fl, eflux)