/
Opaccouls.f
executable file
·121 lines (93 loc) · 4.52 KB
/
Opaccouls.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
c******************************************************************************
c This file contains the functions coulx, coulbf1s, and coulff; these
c routines are used by the continuous opacity routines.
c******************************************************************************
real*8 function coulx (n,freq,z)
c******************************************************************************
c This is function COULX from ATLAS.
c******************************************************************************`
implicit real*8 (a-h,o-z)
real*8 a(6), b(6), c(6)
data a/ 0.9916, 1.105, 1.101, 1.101, 1.102, 1.0986/
data b/ 2.719d13, -2.375d14, -9.863d13, -5.765d13, -3.909d13,
. -2.704d13/
data c/ -2.268d30, 4.077d28, 1.035d28, 4.593d27, 2.371d27,
. 1.229d27/
if (freq .lt. z*z*3.28805d15/dfloat(n*n)) then
coulx = 0.
else
coulx = 2.815d29/freq**3/dfloat(n**5)*z**4
if (n .gt. 6) return
if (n .eq. 1) then
coulx = coulx*coulbf1s(freq,z)
else
coulx = coulx*(a(n)+(b(n)+c(n)*(z*z/freq))*(z*z/freq))
endif
endif
return
end
real*8 function coulbf1s (freq,z)
c******************************************************************************
c This is function COULBF1S from ATLAS.
c******************************************************************************`
implicit real*8 (a-h,o-z)
real*8 gaunt1s(151)
data gaunt1s/
. 0.7973,0.8094,0.8212,0.8328,0.8439,0.8548,0.8653,0.8754,0.8852,
. 0.8946,0.9035,0.9120,0.9201,0.9278,0.9351,0.9420,0.9484,0.9544,
. 0.9601,0.9653,0.9702,0.9745,0.9785,0.9820,0.9852,0.9879,0.9903,
. 0.9922,0.9938,0.9949,0.9957,0.9960,0.9960,0.9957,0.9949,0.9938,
. 0.9923,0.9905,0.9884,0.9859,0.9832,0.9801,0.9767,0.9730,0.9688,
. 0.9645,0.9598,0.9550,0.9499,0.9445,0.9389,0.9330,0.9269,0.9206,
. 0.9140,0.9071,0.9001,0.8930,0.8856,0.8781,0.8705,0.8627,0.8546,
. 0.8464,0.8381,0.8298,0.8213,0.8128,0.8042,0.7954,0.7866,0.7777,
. 0.7685,0.7593,0.7502,0.7410,0.7318,0.7226,0.7134,0.7042,0.6951,
. 0.6859,0.6767,0.6675,0.6584,0.6492,0.6401,0.6310,0.6219,0.6129,
. 0.6039,0.5948,0.5859,0.5769,0.5680,0.5590,0.5502,0.5413,0.5324,
. 0.5236,0.5148,0.5063,0.4979,0.4896,0.4814,0.4733,0.4652,0.4572,
. 0.4493,0.4415,0.4337,0.4261,0.4185,0.4110,0.4035,0.3962,0.3889,
. 0.3818,0.3749,0.3680,0.3611,0.3544,0.3478,0.3413,0.3348,0.3285,
. 0.3222,0.3160,0.3099,0.3039,0.2980,0.2923,0.2866,0.2810,0.2755,
. 0.2701,0.2648,0.2595,0.2544,0.2493,0.2443,0.2394,0.2345,0.2298,
. 0.2251,0.2205,0.2160,0.2115,0.2072,0.2029,0.1987/
if (freq/z**2 .lt. 3.28805d15) then
coulbf1s = 0.
else
elog = log10(freq/z**2/3.28805d15)
i = elog/0.02
i = max(min(i+1,150),1)
coulbf1s = gaunt1s(i) +
. (gaunt1s(i+1)-gaunt1s(i))/0.02*(elog-(i-1)*0.02)
endif
return
end
real*8 function coulff (nz,tlog,freq)
c******************************************************************************
c This is function COULFF from ATLAS.
c******************************************************************************`
implicit real*8 (a-h,o-z)
real*8 z4log(6), a(11,12)
data z4log/ 0., 1.20412, 1.90849, 2.40824, 2.79588, 3.11261/
data a/
. 5.53,5.49,5.46,5.43,5.40,5.25,5.00,4.69,4.48,4.16,3.85,
. 4.91,4.87,4.84,4.80,4.77,4.63,4.40,4.13,3.87,3.52,3.27,
. 4.29,4.25,4.22,4.18,4.15,4.02,3.80,3.57,3.27,2.98,2.70,
. 3.64,3.61,3.59,3.56,3.54,3.41,3.22,2.97,2.70,2.45,2.20,
. 3.00,2.98,2.97,2.95,2.94,2.81,2.65,2.44,2.21,2.01,1.81,
. 2.41,2.41,2.41,2.41,2.41,2.32,2.19,2.02,1.84,1.67,1.50,
. 1.87,1.89,1.91,1.93,1.95,1.90,1.80,1.68,1.52,1.41,1.30,
. 1.33,1.39,1.44,1.49,1.55,1.56,1.51,1.42,1.33,1.25,1.17,
. 0.90,0.95,1.00,1.08,1.17,1.30,1.32,1.30,1.20,1.15,1.11,
. 0.55,0.58,0.62,0.70,0.85,1.01,1.15,1.18,1.15,1.11,1.08,
. 0.33,0.36,0.39,0.46,0.59,0.76,0.97,1.09,1.13,1.10,1.08,
. 0.19,0.21,0.24,0.28,0.38,0.53,0.76,0.96,1.08,1.09,1.09/
gamlog = 10.39638 - tlog/1.15129 + z4log(nz)
igam = max0(min0(int(gamlog+7.),10),1)
hvktlg = (dlog(freq) - tlog)/1.15129 - 20.63764
ihvkt = max0(min0(int(hvktlg+9.),11),1)
p = gamlog - dfloat(igam-7)
q = hvktlg - dfloat(ihvkt-9)
coulff = (1.-p)*((1.-q)*a(igam,ihvkt) + q*a(igam,ihvkt+1)) +
. p*((1.-q)*a(igam+1,ihvkt) + q*a(igam+1,ihvkt+1))
return
end