« back to all changes in this revision
Viewing changes to vendor/SMWidth/hdecay/elw.f
- browse files at revision 276
- compare with another revision
- download diff
- download tarball
- view history from revision 276
- Committer: olivier-mattelaer
- Date: 2018-04-29 08:10:16 UTC
- mfrom: (275.1.80 2.6.2)
- Revision ID: olivier-mattelaer-20180429081016-9nmfvn1er0zjb23o
pass to 2.6.2
- files added:
- Template/NLO/Utilities/open_output_files_dummy_HwU.f
- files removed:
- Template/LO/Source/PDF/Ctq4Fn.f
- vendor/SMWidth
- vendor/SMWidth/hdecay
- vendor/SMWidth/oneloop
- files modified:
- MadSpin/decay.py
added
removed
vendor/SMWidth/hdecay/elw.f
1
*******************************************************************************
2
* *
3
* EWgint *
4
* *
5
* Interpolating the EW grid in gg --> H *
6
* *
7
* November 2009, G. Passarino, C. Sturm, S, Uccirati *
8
* *
9
* *
10
* S.~Actis, G.~Passarino, C.~Sturm and S.~Uccirati, *
11
* Nucl.\ Phys.\ B {\bf 811} (2009) 182 *
12
* [arXiv:0809.3667 [hep-ph]]. *
13
* *
14
* S.~Actis, G.~Passarino, C.~Sturm and S.~Uccirati, *
15
* Phys.\ Lett.\ B {\bf 670} (2008) 12 *
16
* [arXiv:0809.1301 [hep-ph]] *
17
* *
18
* S.~Actis, G.~Passarino, C.~Sturm and S.~Uccirati, *
19
* Phys.\ Lett.\ B {\bf 669} (2008) 62 *
20
* [arXiv:0809.1302 [hep-ph]]. *
21
* *
22
* *
23
* Email: giampiero@to.infn.it *
24
* *
25
*******************************************************************************
26
*
27
double precision function glgl_elw(mtop,mh)
28
*
29
IMPLICIT NONE
30
*
31
INTEGER i,top,gdim
32
REAL*8 u,ui,value
33
REAL*8 vp,vpp,vppp
34
REAL*8 mtop,mh
35
REAL*8 x,x0,x1,x2,y0,y1,y2,a0,a1,a2
36
REAL*8 value0,value1,value2
37
REAL*8 bl(154),cl(154),dl(154)
38
REAL*8 bc(151),cc(151),dc(151)
39
REAL*8 bu(152),cu(152),du(152)
40
*
41
* u value of M_H at which the spline is to be evaluated
42
* top= -1,0,1 lower, central, upper value for m_top
43
*
44
u = mh
45
top = -1
46
gdim= 154
47
CALL FMMsplineSingle(bl,cl,dl,top,gdim)
48
CALL Seval3Single(u,bl,cl,dl,top,gdim,value0,vp,vpp,vppp)
49
top = 0
50
gdim= 151
51
CALL FMMsplineSingle(bc,cc,dc,top,gdim)
52
CALL Seval3Single(u,bc,cc,dc,top,gdim,value1,vp,vpp,vppp)
53
top = 1
54
gdim= 152
55
CALL FMMsplineSingle(bu,cu,du,top,gdim)
56
CALL Seval3Single(u,bu,cu,du,top,gdim,value2,vp,vpp,vppp)
57
x = mtop
58
x0 = 171.06d0
59
x1 = 172.64d0
60
x2 = 174.22d0
61
y0 = value0
62
y1 = value1
63
y2 = value2
64
A0=(X-X1)*(X-X2)/(X0-X1)/(X0-X2)
65
A1=(X-X0)*(X-X2)/(X1-X0)/(X1-X2)
66
A2=(X-X0)*(X-X1)/(X2-X0)/(X2-X1)
67
glgl_elw=(A0*Y0+A1*Y1+A2*Y2)/100.d0
68
RETURN
69
END
70
*
71
*-----------------------------------------------------------------------
72
*
73
c CONTAINS
74
*
75
SUBROUTINE FMMsplineSingle(b,c,d,top,gdim)
76
*
77
* ---------------------------------------------------------------------------
78
*
79
* PURPOSE - Compute the coefficients b,c,d for a cubic interpolating spline
80
* so that the interpolated value is given by
81
* s(x) = y(k) + b(k)*(x-x(k)) + c(k)*(x-x(k))**2 + d(k)*(x-x(k))**3
82
* when x(k) <= x <= x(k+1)
83
* The end conditions match the third derivatives of the interpolated curve to
84
* the third derivatives of the unique polynomials thru the first four and
85
* last four points.
86
* Use Seval or Seval3 to evaluate the spline.
87
*
88
INTEGER k,n,i,top,gdim,l
89
*
90
REAL*8 xl(154),yl(154)
91
REAL*8 xc(151),yc(151)
92
REAL*8 xu(152),yu(152)
93
REAL*8 x(154),y(154)
94
*
95
REAL*8 b(gdim)
96
* linear coeff
97
*
98
REAL*8 c(gdim)
99
* quadratic coeff.
100
*
101
REAL*8 d(gdim)
102
* cubic coeff.
103
*
104
REAL*8 t
105
REAL*8 ZERO, TWO, THREE
106
PARAMETER(ZERO=0.0, TWO=2.0, THREE=3.0)
107
*
108
* The grid
109
*
110
DATA (xl(i),i=1,154)/100.0d0,110.0d0,120.0d0,130.0d0,
111
# 140.0d0,145.0d0,150.0d0,151.0d0,
112
# 152.0d0,153.0d0,154.0d0,155.0d0,
113
# 156.0d0,157.0d0,158.0d0,159.0d0,
114
# 160.0d0,161.0d0,162.0d0,163.0d0,
115
# 164.0d0,165.0d0,166.0d0,
116
# 167.0d0,168.0d0,169.0d0,170.0d0,
117
# 171.0d0,172.0d0,173.0d0,174.0d0,
118
# 175.0d0,176.0d0,177.0d0,178.0d0,
119
# 179.0d0,180.0d0,181.0d0,182.0d0,
120
# 183.0d0,184.0d0,185.0d0,186.0d0,
121
# 187.0d0,188.0d0,189.0d0,190.0d0,
122
# 195.0d0,197.0d0,200.0d0,210.0d0,
123
# 220.0d0,230.0d0,240.0d0,250.0d0,
124
# 260.0d0,270.0d0,280.0d0,290.0d0,
125
# 300.0d0,310.0d0,320.0d0,330.0d0,
126
# 335.0d0,340.0d0,341.0d0,342.0d0,
127
# 343.0d0,344.0d0,345.0d0,345.3d0,
128
# 345.5d0,345.8d0,346.0d0,347.0d0,
129
# 348.0d0,349.0d0,349.5d0,349.75d0,
130
# 350.0d0,351.0d0,352.0d0,353.0d0,
131
# 354.0d0,355.0d0,356.0d0,357.0d0,
132
# 358.0d0,359.0d0,360.0d0,370.0d0,
133
# 380.0d0,390.0d0,400.0d0,410.0d0,
134
# 420.0d0,430.0d0,440.0d0,450.0d0,
135
# 460.0d0,470.0d0,480.0d0,490.0d0,
136
# 500.0d0,510.0d0,520.0d0,530.0d0,
137
# 540.0d0,550.0d0,560.0d0,570.0d0,
138
# 580.0d0,590.0d0,600.0d0,610.0d0,
139
# 620.0d0,630.0d0,640.0d0,650.0d0,
140
# 660.0d0,670.0d0,680.0d0,690.0d0,
141
# 700.0d0,710.0d0,720.0d0,730.0d0,
142
# 740.0d0,750.0d0,760.0d0,770.0d0,
143
# 780.0d0,790.0d0,800.0d0,810.0d0,
144
# 820.0d0,830.0d0,840.0d0,850.0d0,
145
# 860.0d0,870.0d0,880.0d0,890.0d0,
146
# 900.0d0,910.0d0,920.0d0,930.0d0,
147
# 940.0d0,950.0d0,960.0d0,970.0d0,
148
# 980.0d0,990.0d0,1000.0d0/
149
*
150
DATA (yl(i),i=1,154)/4.179467154d0,4.542088751d0,
151
# 4.938014960d0,5.335142457d0,
152
# 5.707514034d0,5.868498547d0,5.978310949d0,5.986920781d0,
153
# 5.987460634d0,5.977179301d0,5.951667414d0,5.904562092d0,
154
# 5.826081147d0,5.700679828d0,5.506218392d0,5.220567269d0,
155
# 4.861752146d0,4.527674442d0,4.222683291d0,3.880374785d0,
156
# 3.529354583d0,3.204116062d0,2.915522583d0,
157
# 2.662502577d0,2.440241875d0,2.243618198d0,2.067597988d0,
158
# 1.907471856d0,1.759641285d0,1.620735593d0,1.487565904d0,
159
# 1.355817302d0,1.222813267d0,1.081630216d0,0.926313642d0,
160
# 0.747981132d0,0.538696787d0,0.306541839d0,0.085775708d0,
161
# -0.092817756d0,-0.259991178d0,-0.443225249d0,-0.632350007d0,
162
# -0.811227830d0,-0.975263093d0,-1.126478218d0,-1.262721067d0,
163
# -1.757119775d0,-1.895339688d0,-2.057961996d0,-2.370870160d0,
164
# -2.488446977d0,-2.490894901d0,-2.433191749d0,-2.359142379d0,
165
# -2.260728908d0,-2.151155918d0,-2.036327859d0,-1.933794725d0,
166
# -1.862582453d0,-1.802070746d0,-1.812386184d0,-1.938240382d0,
167
# -2.108006605d0,-2.488306576d0,-2.655000325d0,-3.016952793d0,
168
# -3.680237527d0,-3.870995189d0,-3.959387796d0,-3.943637421d0,
169
# -3.982342165d0,-3.967290244d0,-4.002598626d0,-4.008430095d0,
170
# -4.010824003d0,-4.021419122d0,-4.018822218d0,-4.040141549d0,
171
# -4.024343762d0,-4.015666188d0,-3.965077563d0,-3.968563552d0,
172
# -3.963667776d0,-3.942755435d0,-3.932719674d0,-3.874191897d0,
173
# -3.876081328d0,-3.833375576d0,-3.796241844d0,-3.406242887d0,
174
# -3.088048772d0,-2.738617943d0,-2.446609831d0,-2.131731318d0,
175
# -1.878917579d0,-1.589864293d0,-1.358374992d0,-1.103676110d0,
176
# -0.884448436d0,-0.679633948d0,-0.472673942d0,-0.281736254d0,
177
# -0.108375079d0,0.093529154d0,0.257607461d0,0.430152846d0,
178
# 0.578396867d0,0.763320409d0,0.934522954d0,1.075378845d0,
179
# 1.249326065d0,1.426785290d0,1.535796936d0,1.715981764d0,
180
# 1.871449481d0,2.020126942d0,2.154595978d0,2.322995662d0,
181
# 2.456079279d0,2.608032219d0,2.749902433d0,2.905996823d0,
182
# 3.064612486d0,3.205141310d0,3.344347828d0,3.506879165d0,
183
# 3.660117972d0,3.783439131d0,3.947284161d0,4.097179155d0,
184
# 4.248755764d0,4.406763712d0,4.563193913d0,4.705927695d0,
185
# 4.853626937d0,5.015246327d0,5.167579841d0,5.341215028d0,
186
# 5.512112004d0,5.650306940d0,5.805087521d0,5.980800439d0,
187
# 6.125262600d0,6.272280377d0,6.477594572d0,6.625428154d0,
188
# 6.810107710d0,6.954753946d0,7.131768272d0,7.304668734d0,
189
# 7.458544663d0,7.653458644d0,7.821154640d0/
190
*
191
DATA (xc(i),i=1,151)/100.0d0,110.0d0,120.0d0,130.0d0,140.0d0,
192
# 145.0d0,150.0d0,151.0d0,152.0d0,153.0d0,154.0d0,155.0d0,
193
# 156.0d0,157.0d0,158.0d0,159.0d0,160.0d0,161.0d0,162.0d0,
194
# 163.0d0,164.0d0,165.0d0,166.0d0,167.0d0,168.0d0,169.0d0,
195
# 170.0d0,171.0d0,172.0d0,173.0d0,174.0d0,175.0d0,176.0d0,
196
# 177.0d0,178.0d0,179.0d0,180.0d0,181.0d0,182.0d0,183.0d0,
197
# 184.0d0,185.0d0,186.0d0,187.0d0,188.0d0,
198
# 189.0d0,190.0d0,195.0d0,197.0d0,200.0d0,
199
# 210.0d0,220.0d0,230.0d0,240.0d0,250.0d0,
200
# 260.0d0,270.0d0,280.0d0,290.0d0,300.0d0,
201
# 320.0d0,330.0d0,335.0d0,340.0d0,341.0d0,
202
# 342.0d0,343.0d0,344.0d0,345.0d0,345.3d0,
203
# 345.5d0,345.8d0,346.0d0,347.0d0,348.0d0,
204
# 349.0d0,350.0d0,351.0d0,352.0d0,353.0d0,
205
# 354.0d0,355.0d0,356.0d0,357.0d0,358.0d0,
206
# 359.0d0,360.0d0,370.0d0,380.0d0,390.0d0,
207
# 400.0d0,410.0d0,420.0d0,430.0d0,440.0d0,
208
# 450.0d0,460.0d0,470.0d0,480.0d0,490.0d0,
209
# 500.0d0,510.0d0,520.0d0,530.0d0,540.0d0,
210
# 550.0d0,560.0d0,570.0d0,580.0d0,590.0d0,
211
# 600.0d0,610.0d0,620.0d0,630.0d0,640.0d0,
212
# 650.0d0,660.0d0,670.0d0,680.0d0,690.0d0,
213
# 700.0d0,710.0d0,720.0d0,730.0d0,740.0d0,
214
# 750.0d0,760.0d0,770.0d0,780.0d0,790.0d0,
215
# 800.0d0,810.0d0,820.0d0,830.0d0,840.0d0,
216
# 850.0d0,860.0d0,870.0d0,880.0d0,890.0d0,
217
# 900.0d0,910.0d0,920.0d0,930.0d0,940.0d0,
218
# 950.0d0,960.0d0,970.0d0,980.0d0,990.0d0, 1000.0d0/
219
*
220
DATA (yc(i),i=1,151)/4.183581334d0,4.545978356d0,
221
# 4.941666452d0,5.338524432d0,
222
# 5.710552332d0,5.871305639d0,5.980798334d0,5.989325219d0,
223
# 5.989771990d0,5.979384655d0,5.953749466d0,5.906497163d0,
224
# 5.827835443d0,5.702202128d0,5.507427195d0,5.221327470d0,
225
# 4.861849085d0,4.526867019d0,4.220897196d0,3.877727401d0,
226
# 3.525995801d0,3.200170282d0,2.911100907d0,2.657698720d0,
227
# 2.435102574d0,2.238185021d0,2.061849416d0,1.901483110d0,
228
# 1.753411931d0,1.614301638d0,1.480965572d0,1.348997661d0,
229
# 1.215803604d0,1.074412250d0,0.919019083d0,0.740332966d0,
230
# 0.530749433d0,0.298400521d0,0.077023976d0,-0.102260900d0,
231
# -0.270062463d0,-0.453870620d0,-0.643469731d0,-0.822931558d0,
232
# -0.987443985d0,-1.138691526d0,-1.275416891d0,-1.770985289d0,
233
# -1.910067996d0,-2.073780511d0,-2.387106650d0,-2.506786934d0,
234
# -2.507793900d0,-2.451529121d0,-2.379171774d0,-2.280179838d0,
235
# -2.173130895d0,-2.052273974d0,-1.953341364d0,-1.870193084d0,
236
# -1.795179493d0,-1.874006942d0,-1.991907667d0,-2.212751797d0,
237
# -2.281770963d0,-2.366744787d0,-2.479318144d0,-2.637428316d0,
238
# -2.936192030d0,-3.392136424d0,-3.518873910d0,-3.682965171d0,
239
# -3.718363034d0,-3.872185129d0,-3.990484562d0,-4.003435184d0,
240
# -4.047831469d0,-4.056871742d0,-4.096278566d0,-4.096147175d0,
241
# -4.069660494d0,-4.092229117d0,-4.050384852d0,-4.043310923d0,
242
# -4.021791606d0,-4.007601020d0,-3.922141199d0,-3.576132878d0,
243
# -3.231059452d0,-2.868110119d0,-2.588448565d0,-2.286823317d0,
244
# -1.981203740d0,-1.697922238d0,-1.465949078d0,-1.221405908d0,
245
# -0.997896019d0,-0.787814864d0,-0.603630999d0,-0.385428569d0,
246
# -0.191489219d0,0.012545468d0,0.180344393d0,0.350524042d0,
247
# 0.535367220d0,0.686993292d0,0.894786143d0,1.012194233d0,
248
# 1.192054257d0,1.364029746d0,1.515391057d0,1.655859873d0,
249
# 1.814367096d0,1.961412767d0,2.115811847d0,2.260602653d0,
250
# 2.429881738d0,2.579117247d0,2.705063788d0,2.852512294d0,
251
# 3.013412959d0,3.190250179d0,3.334985916d0,3.458011579d0,
252
# 3.596519369d0,3.772420987d0,3.927207591d0,4.083443432d0,
253
# 4.220565084d0,4.367619511d0,4.540226035d0,4.671279216d0,
254
# 4.822978517d0,5.004941301d0,5.133890543d0,5.317216628d0,
255
# 5.516475363d0,5.649275363d0,5.771911413d0,5.988071980d0,
256
# 6.100250110d0,6.265613527d0,6.437260016d0,6.622288364d0,
257
# 6.781564747d0,6.953709972d0,7.135417771d0,7.291456432d0,
258
# 7.475317109d0,7.623304078d0,7.829501710d0/
259
*
260
DATA (xu(i),i=1,152)/100.0d0,110.0d0,120.0d0,130.0d0,
261
# 140.0d0,145.0d0,150.0d0,151.0d0,
262
# 152.0d0,153.0d0,154.0d0,155.0d0,
263
# 156.0d0,157.0d0,158.0d0,159.0d0,
264
# 160.0d0,161.0d0,162.0d0,163.0d0,
265
# 164.0d0,165.0d0,166.0d0,167.0d0,
266
# 168.0d0,169.0d0,170.0d0,171.0d0,
267
# 172.0d0,173.0d0,174.0d0,175.0d0,
268
# 176.0d0,177.0d0,178.0d0,179.0d0,
269
# 180.0d0,181.0d0,182.0d0,183.0d0,
270
# 184.0d0,185.0d0,186.0d0,187.0d0,
271
# 188.0d0,189.0d0,190.0d0,195.0d0,
272
# 197.0d0,200.0d0,210.0d0,220.0d0,
273
# 230.0d0,240.0d0,250.0d0,260.0d0,
274
# 270.0d0,280.0d0,290.0d0,300.0d0,
275
# 310.0d0,320.0d0,330.0d0,335.0d0,
276
# 340.0d0,341.0d0,342.0d0,343.0d0,
277
# 344.0d0,345.0d0,345.3d0,345.5d0,
278
# 345.8d0,346.0d0,347.0d0,348.0d0,
279
# 349.0d0,350.0d0,351.0d0,352.0d0,
280
# 353.0d0,354.0d0,355.0d0,356.0d0,
281
# 357.0d0,358.0d0,359.0d0,360.0d0,
282
# 370.0d0,380.0d0,390.0d0,400.0d0,
283
# 410.0d0,420.0d0,430.0d0,440.0d0,
284
# 450.0d0,460.0d0,470.0d0,480.0d0,
285
# 490.0d0,500.0d0,510.0d0,520.0d0,
286
# 530.0d0,540.0d0,550.0d0,560.0d0,
287
# 570.0d0,580.0d0,590.0d0,600.0d0,
288
# 610.0d0,620.0d0,630.0d0,640.0d0,
289
# 650.0d0,660.0d0,670.0d0,680.0d0,
290
# 690.0d0,700.0d0,710.0d0,720.0d0,
291
# 730.0d0,740.0d0,750.0d0,760.0d0,
292
# 770.0d0,780.0d0,790.0d0,800.0d0,
293
# 810.0d0,820.0d0,830.0d0,840.0d0,
294
# 850.0d0,860.0d0,870.0d0,880.0d0,
295
# 890.0d0,900.0d0,910.0d0,920.0d0,
296
# 930.0d0,940.0d0,950.0d0,960.0d0,
297
# 970.0d0,980.0d0,990.0d0,1000.0d0/
298
*
299
DATA (yu(i),i=1,152)/4.187720723d0,4.549885500d0,
300
# 4.945324580d0,5.341899188d0,
301
# 5.713567160d0,5.874080782d0,5.983246138d0,5.991688858d0,
302
# 5.992041524d0,5.981547445d0,5.955788600d0,5.908389532d0,
303
# 5.829548098d0,5.703685183d0,5.508601342d0,5.222061322d0,
304
# 4.861934021d0,4.526068905d0,4.219144123d0,3.875133708d0,
305
# 3.522706643d0,3.196306140d0,2.906767625d0,2.652991136d0,
306
# 2.430066986d0,2.232861152d0,2.056219593d0,1.895596836d0,
307
# 1.747309568d0,1.607964744d0,1.474493488d0,1.342295851d0,
308
# 1.208919594d0,1.067309617d0,0.911812113d0,0.732879892d0,
309
# 0.522870581d0,0.290366470d0,0.068493581d0,-0.111577414d0,
310
# -0.279962552d0,-0.464332982d0,-0.654447961d0,-0.834336934d0,
311
# -0.999445210d0,-1.150871949d0,-1.287684579d0,-1.784688615d0,
312
# -1.924511090d0,-2.089081933d0,-2.402982024d0,-2.524611659d0,
313
# -2.525387652d0,-2.470417238d0,-2.399507117d0,-2.298953258d0,
314
# -2.192459462d0,-2.069700408d0,-1.973189437d0,-1.880904758d0,
315
# -1.807288585d0,-1.785086586d0,-1.826974607d0,-1.915990447d0,
316
# -2.062701208d0,-2.104793465d0,-2.152577512d0,-2.208540347d0,
317
# -2.276262386d0,-2.361739355d0,-2.391880716d0,-2.413272954d0,
318
# -2.447493257d0,-2.471820336d0,-2.622996766d0,-2.885017381d0,
319
# -3.696189261d0,-3.893536695d0,-4.011725974d0,-4.074946087d0,
320
# -4.137470131d0,-4.089349430d0,-4.136099918d0,-4.123934918d0,
321
# -4.136446777d0,-4.132255956d0,-4.118820717d0,-4.088009318d0,
322
# -3.788430916d0,-3.406565873d0,-3.050120268d0,-2.724817324d0,
323
# -2.410462214d0,-2.120474142d0,-1.845522666d0,-1.583137592d0,
324
# -1.337585290d0,-1.094144564d0,-0.897029424d0,-0.671551439d0,
325
# -0.478508120d0,-0.269184704d0,-0.072526771d0,0.115732816d0,
326
# 0.274371119d0,0.457987609d0,0.619147501d0,0.809852863d0,
327
# 0.968131131d0,1.121473977d0,1.290603896d0,1.462164252d0,
328
# 1.583836738d0,1.785815006d0,1.912473989d0,2.068564737d0,
329
# 2.235729083d0,2.398009567d0,2.538161792d0,2.681468958d0,
330
# 2.803535096d0,2.966400538d0,3.143956306d0,3.310124390d0,
331
# 3.429746138d0,3.580977910d0,3.727340256d0,3.890426351d0,
332
# 4.045983092d0,4.203320822d0,4.348209516d0,4.497115304d0,
333
# 4.654911988d0,4.823398533d0,4.983830013d0,5.127442851d0,
334
# 5.265056876d0,5.448819268d0,5.642707192d0,5.785441771d0,
335
# 5.936113103d0,6.108286154d0,6.258324908d0,6.415575321d0,
336
# 6.591607190d0,6.761124005d0,6.951304838d0,7.117437285d0,
337
# 7.283053264d0,7.453490877d0,7.637739897d0,7.799659928d0/
338
*
339
IF(top.eq.-1) THEN
340
n= 154
341
DO l=1,154
342
x(l)= xl(l)
343
y(l)= yl(l)
344
ENDDO
345
ELSEIF(top.eq.0) THEN
346
n= 151
347
DO l=1,151
348
x(l)= xc(l)
349
y(l)= yc(l)
350
ENDDO
351
ELSEIF(top.eq.1) THEN
352
n= 152
353
DO l=1,152
354
x(l)= xu(l)
355
y(l)= yu(l)
356
ENDDO
357
ENDIF
358
359
*.....Set up tridiagonal system.........................................
360
* b=diagonal, d=offdiagonal, c=right-hand side
361
*
362
d(1)= x(2)-x(1)
363
c(2)= (y(2)-y(1))/d(1)
364
DO k= 2,n-1
365
d(k)= x(k+1)-x(k)
366
b(k)= TWO*(d(k-1)+d(k))
367
c(k+1)= (y(k+1)-y(k))/d(k)
368
c(k)= c(k+1)-c(k)
369
END DO
370
*
371
*.....End conditions. third derivatives at x(1) and x(n) obtained
372
* from divided differences.......................................
373
*
374
b(1)= -d(1)
375
b(n)= -d(n-1)
376
c(1)= ZERO
377
c(n)= ZERO
378
IF (n.gt.3) THEN
379
c(1)= c(3)/(x(4)-x(2))-c(2)/(x(3)-x(1))
380
c(n)= c(n-1)/(x(n)-x(n-2))-c(n-2)/(x(n-1)-x(n-3))
381
c(1)= c(1)*d(1)*d(1)/(x(4)-x(1))
382
c(n)= -c(n)*d(n-1)*d(n-1)/(x(n)-x(n-3))
383
END IF
384
*
385
DO k=2,n ! forward elimination
386
t= d(k-1)/b(k-1)
387
b(k)= b(k)-t*d(k-1)
388
c(k)= c(k)-t*c(k-1)
389
END DO
390
*
391
c(n)= c(n)/b(n)
392
*
393
* back substitution ( makes c the sigma of text)
394
*
395
DO k=n-1,1,-1
396
c(k)= (c(k)-d(k)*c(k+1))/b(k)
397
END DO
398
*
399
*.....Compute polynomial coefficients...................................
400
*
401
b(n)= (y(n)-y(n-1))/d(n-1)+d(n-1)*(c(n-1)+c(n)+c(n))
402
DO k=1,n-1
403
b(k)= (y(k+1)-y(k))/d(k)-d(k)*(c(k+1)+c(k)+c(k))
404
d(k)= (c(k+1)-c(k))/d(k)
405
c(k)= THREE*c(k)
406
END DO
407
c(n)= THREE*c(n)
408
d(n)= d(n-1)
409
*
410
RETURN
411
*
412
END
413
*
414
*------------------------------------------------------------------------
415
*
416
SUBROUTINE Seval3Single(u,b,c,d,top,gdim,f,fp,fpp,fppp)
417
*
418
* ---------------------------------------------------------------------------
419
*
420
* PURPOSE - Evaluate the cubic spline function
421
* Seval=y(i)+b(i)*(u-x(i))+c(i)*(u-x(i))**2+d(i)*(u-x(i))**3
422
* where x(i) <= u < x(i+1)
423
*
424
* NOTES- if u<x(1), i=1 is used;if u>x(n), i=n is used
425
*
426
REAL*8 u
427
* abscissa at which the spline is to be evaluated
428
*
429
INTEGER j,k,n,l,top,gdim
430
*
431
REAL*8 xl(154),yl(154)
432
REAL*8 xc(151),yc(151)
433
REAL*8 xu(152),yu(152)
434
REAL*8 x(154),y(154)
435
REAL*8 b(gdim),c(gdim),d(gdim)
436
* linear,quadratic,cubic coeff
437
*
438
REAL*8 f,fp,fpp,fppp
439
* function, 1st,2nd,3rd deriv
440
*
441
INTEGER i
442
REAL*8 dx
443
REAL*8 ZERO, TWO, THREE
444
PARAMETER(ZERO=0.0, TWO=2.0, THREE=3.0)
445
SAVE i
446
DATA i/1/
447
*
448
* The grid
449
*
450
DATA (xl(l),l= 1,154)/100.0d0,110.0d0,120.0d0,130.0d0,
451
# 140.0d0,145.0d0,150.0d0,151.0d0,
452
# 152.0d0,153.0d0,154.0d0,155.0d0,
453
# 156.0d0,157.0d0,158.0d0,159.0d0,
454
# 160.0d0,161.0d0,162.0d0,163.0d0,
455
# 164.0d0,165.0d0,166.0d0,
456
# 167.0d0,168.0d0,169.0d0,170.0d0,
457
# 171.0d0,172.0d0,173.0d0,174.0d0,
458
# 175.0d0,176.0d0,177.0d0,178.0d0,
459
# 179.0d0,180.0d0,181.0d0,182.0d0,
460
# 183.0d0,184.0d0,185.0d0,186.0d0,
461
# 187.0d0,188.0d0,189.0d0,190.0d0,
462
# 195.0d0,197.0d0,200.0d0,210.0d0,
463
# 220.0d0,230.0d0,240.0d0,250.0d0,
464
# 260.0d0,270.0d0,280.0d0,290.0d0,
465
# 300.0d0,310.0d0,320.0d0,330.0d0,
466
# 335.0d0,340.0d0,341.0d0,342.0d0,
467
# 343.0d0,344.0d0,345.0d0,345.3d0,
468
# 345.5d0,345.8d0,346.0d0,347.0d0,
469
# 348.0d0,349.0d0,349.5d0,349.75d0,
470
# 350.0d0,351.0d0,352.0d0,353.0d0,
471
# 354.0d0,355.0d0,356.0d0,357.0d0,
472
# 358.0d0,359.0d0,360.0d0,370.0d0,
473
# 380.0d0,390.0d0,400.0d0,410.0d0,
474
# 420.0d0,430.0d0,440.0d0,450.0d0,
475
# 460.0d0,470.0d0,480.0d0,490.0d0,
476
# 500.0d0,510.0d0,520.0d0,530.0d0,
477
# 540.0d0,550.0d0,560.0d0,570.0d0,
478
# 580.0d0,590.0d0,600.0d0,610.0d0,
479
# 620.0d0,630.0d0,640.0d0,650.0d0,
480
# 660.0d0,670.0d0,680.0d0,690.0d0,
481
# 700.0d0,710.0d0,720.0d0,730.0d0,
482
# 740.0d0,750.0d0,760.0d0,770.0d0,
483
# 780.0d0,790.0d0,800.0d0,810.0d0,
484
# 820.0d0,830.0d0,840.0d0,850.0d0,
485
# 860.0d0,870.0d0,880.0d0,890.0d0,
486
# 900.0d0,910.0d0,920.0d0,930.0d0,
487
# 940.0d0,950.0d0,960.0d0,970.0d0,
488
# 980.0d0,990.0d0,1000.0d0/
489
*
490
DATA (yl(l),l=1,154)/4.179467154d0,4.542088751d0,
491
# 4.938014960d0,5.335142457d0,
492
# 5.707514034d0,5.868498547d0,5.978310949d0,5.986920781d0,
493
# 5.987460634d0,5.977179301d0,5.951667414d0,5.904562092d0,
494
# 5.826081147d0,5.700679828d0,5.506218392d0,5.220567269d0,
495
# 4.861752146d0,4.527674442d0,4.222683291d0,3.880374785d0,
496
# 3.529354583d0,3.204116062d0,2.915522583d0,
497
# 2.662502577d0,2.440241875d0,2.243618198d0,2.067597988d0,
498
# 1.907471856d0,1.759641285d0,1.620735593d0,1.487565904d0,
499
# 1.355817302d0,1.222813267d0,1.081630216d0,0.926313642d0,
500
# 0.747981132d0,0.538696787d0,0.306541839d0,0.085775708d0,
501
# -0.092817756d0,-0.259991178d0,-0.443225249d0,-0.632350007d0,
502
# -0.811227830d0,-0.975263093d0,-1.126478218d0,-1.262721067d0,
503
# -1.757119775d0,-1.895339688d0,-2.057961996d0,-2.370870160d0,
504
# -2.488446977d0,-2.490894901d0,-2.433191749d0,-2.359142379d0,
505
# -2.260728908d0,-2.151155918d0,-2.036327859d0,-1.933794725d0,
506
# -1.862582453d0,-1.802070746d0,-1.812386184d0,-1.938240382d0,
507
# -2.108006605d0,-2.488306576d0,-2.655000325d0,-3.016952793d0,
508
# -3.680237527d0,-3.870995189d0,-3.959387796d0,-3.943637421d0,
509
# -3.982342165d0,-3.967290244d0,-4.002598626d0,-4.008430095d0,
510
# -4.010824003d0,-4.021419122d0,-4.018822218d0,-4.040141549d0,
511
# -4.024343762d0,-4.015666188d0,-3.965077563d0,-3.968563552d0,
512
# -3.963667776d0,-3.942755435d0,-3.932719674d0,-3.874191897d0,
513
# -3.876081328d0,-3.833375576d0,-3.796241844d0,-3.406242887d0,
514
# -3.088048772d0,-2.738617943d0,-2.446609831d0,-2.131731318d0,
515
# -1.878917579d0,-1.589864293d0,-1.358374992d0,-1.103676110d0,
516
# -0.884448436d0,-0.679633948d0,-0.472673942d0,-0.281736254d0,
517
# -0.108375079d0,0.093529154d0,0.257607461d0,0.430152846d0,
518
# 0.578396867d0,0.763320409d0,0.934522954d0,1.075378845d0,
519
# 1.249326065d0,1.426785290d0,1.535796936d0,1.715981764d0,
520
# 1.871449481d0,2.020126942d0,2.154595978d0,2.322995662d0,
521
# 2.456079279d0,2.608032219d0,2.749902433d0,2.905996823d0,
522
# 3.064612486d0,3.205141310d0,3.344347828d0,3.506879165d0,
523
# 3.660117972d0,3.783439131d0,3.947284161d0,4.097179155d0,
524
# 4.248755764d0,4.406763712d0,4.563193913d0,4.705927695d0,
525
# 4.853626937d0,5.015246327d0,5.167579841d0,5.341215028d0,
526
# 5.512112004d0,5.650306940d0,5.805087521d0,5.980800439d0,
527
# 6.125262600d0,6.272280377d0,6.477594572d0,6.625428154d0,
528
# 6.810107710d0,6.954753946d0,7.131768272d0,7.304668734d0,
529
# 7.458544663d0,7.653458644d0,7.821154640d0/
530
*
531
DATA (xc(l),l= 1,151)/100.0d0,110.0d0,120.0d0,130.0d0,140.0d0,
532
# 145.0d0,150.0d0,151.0d0,152.0d0,153.0d0,154.0d0,155.0d0,
533
# 156.0d0,157.0d0,158.0d0,159.0d0,160.0d0,161.0d0,162.0d0,
534
# 163.0d0,164.0d0,165.0d0,166.0d0,167.0d0,168.0d0,169.0d0,
535
# 170.0d0,171.0d0,172.0d0,173.0d0,174.0d0,175.0d0,176.0d0,
536
# 177.0d0,178.0d0,179.0d0,180.0d0,181.0d0,182.0d0,183.0d0,
537
# 184.0d0,185.0d0,186.0d0,187.0d0,188.0d0,
538
# 189.0d0,190.0d0,195.0d0,197.0d0,200.0d0,
539
# 210.0d0,220.0d0,230.0d0,240.0d0,250.0d0,
540
# 260.0d0,270.0d0,280.0d0,290.0d0,300.0d0,
541
# 320.0d0,330.0d0,335.0d0,340.0d0,341.0d0,
542
# 342.0d0,343.0d0,344.0d0,345.0d0,345.3d0,
543
# 345.5d0,345.8d0,346.0d0,347.0d0,348.0d0,
544
# 349.0d0,350.0d0,351.0d0,352.0d0,353.0d0,
545
# 354.0d0,355.0d0,356.0d0,357.0d0,358.0d0,
546
# 359.0d0,360.0d0,370.0d0,380.0d0,390.0d0,
547
# 400.0d0,410.0d0,420.0d0,430.0d0,440.0d0,
548
# 450.0d0,460.0d0,470.0d0,480.0d0,490.0d0,
549
# 500.0d0,510.0d0,520.0d0,530.0d0,540.0d0,
550
# 550.0d0,560.0d0,570.0d0,580.0d0,590.0d0,
551
# 600.0d0,610.0d0,620.0d0,630.0d0,640.0d0,
552
# 650.0d0,660.0d0,670.0d0,680.0d0,690.0d0,
553
# 700.0d0,710.0d0,720.0d0,730.0d0,740.0d0,
554
# 750.0d0,760.0d0,770.0d0,780.0d0,790.0d0,
555
# 800.0d0,810.0d0,820.0d0,830.0d0,840.0d0,
556
# 850.0d0,860.0d0,870.0d0,880.0d0,890.0d0,
557
# 900.0d0,910.0d0,920.0d0,930.0d0,940.0d0,
558
# 950.0d0,960.0d0,970.0d0,980.0d0,990.0d0, 1000.0d0/
559
*
560
DATA (yc(l),l= 1,151)/4.183581334d0,4.545978356d0,
561
# 4.941666452d0,5.338524432d0,
562
# 5.710552332d0,5.871305639d0,5.980798334d0,5.989325219d0,
563
# 5.989771990d0,5.979384655d0,5.953749466d0,5.906497163d0,
564
# 5.827835443d0,5.702202128d0,5.507427195d0,5.221327470d0,
565
# 4.861849085d0,4.526867019d0,4.220897196d0,3.877727401d0,
566
# 3.525995801d0,3.200170282d0,2.911100907d0,2.657698720d0,
567
# 2.435102574d0,2.238185021d0,2.061849416d0,1.901483110d0,
568
# 1.753411931d0,1.614301638d0,1.480965572d0,1.348997661d0,
569
# 1.215803604d0,1.074412250d0,0.919019083d0,0.740332966d0,
570
# 0.530749433d0,0.298400521d0,0.077023976d0,-0.102260900d0,
571
# -0.270062463d0,-0.453870620d0,-0.643469731d0,-0.822931558d0,
572
# -0.987443985d0,-1.138691526d0,-1.275416891d0,-1.770985289d0,
573
# -1.910067996d0,-2.073780511d0,-2.387106650d0,-2.506786934d0,
574
# -2.507793900d0,-2.451529121d0,-2.379171774d0,-2.280179838d0,
575
# -2.173130895d0,-2.052273974d0,-1.953341364d0,-1.870193084d0,
576
# -1.795179493d0,-1.874006942d0,-1.991907667d0,-2.212751797d0,
577
# -2.281770963d0,-2.366744787d0,-2.479318144d0,-2.637428316d0,
578
# -2.936192030d0,-3.392136424d0,-3.518873910d0,-3.682965171d0,
579
# -3.718363034d0,-3.872185129d0,-3.990484562d0,-4.003435184d0,
580
# -4.047831469d0,-4.056871742d0,-4.096278566d0,-4.096147175d0,
581
# -4.069660494d0,-4.092229117d0,-4.050384852d0,-4.043310923d0,
582
# -4.021791606d0,-4.007601020d0,-3.922141199d0,-3.576132878d0,
583
# -3.231059452d0,-2.868110119d0,-2.588448565d0,-2.286823317d0,
584
# -1.981203740d0,-1.697922238d0,-1.465949078d0,-1.221405908d0,
585
# -0.997896019d0,-0.787814864d0,-0.603630999d0,-0.385428569d0,
586
# -0.191489219d0,0.012545468d0,0.180344393d0,0.350524042d0,
587
# 0.535367220d0,0.686993292d0,0.894786143d0,1.012194233d0,
588
# 1.192054257d0,1.364029746d0,1.515391057d0,1.655859873d0,
589
# 1.814367096d0,1.961412767d0,2.115811847d0,2.260602653d0,
590
# 2.429881738d0,2.579117247d0,2.705063788d0,2.852512294d0,
591
# 3.013412959d0,3.190250179d0,3.334985916d0,3.458011579d0,
592
# 3.596519369d0,3.772420987d0,3.927207591d0,4.083443432d0,
593
# 4.220565084d0,4.367619511d0,4.540226035d0,4.671279216d0,
594
# 4.822978517d0,5.004941301d0,5.133890543d0,5.317216628d0,
595
# 5.516475363d0,5.649275363d0,5.771911413d0,5.988071980d0,
596
# 6.100250110d0,6.265613527d0,6.437260016d0,6.622288364d0,
597
# 6.781564747d0,6.953709972d0,7.135417771d0,7.291456432d0,
598
# 7.475317109d0,7.623304078d0,7.829501710d0/
599
*
600
DATA (xu(l),l=1,152)/100.0d0,110.0d0,120.0d0,130.0d0,
601
# 140.0d0,145.0d0,150.0d0,151.0d0,
602
# 152.0d0,153.0d0,154.0d0,155.0d0,
603
# 156.0d0,157.0d0,158.0d0,159.0d0,
604
# 160.0d0,161.0d0,162.0d0,163.0d0,
605
# 164.0d0,165.0d0,166.0d0,167.0d0,
606
# 168.0d0,169.0d0,170.0d0,171.0d0,
607
# 172.0d0,173.0d0,174.0d0,175.0d0,
608
# 176.0d0,177.0d0,178.0d0,179.0d0,
609
# 180.0d0,181.0d0,182.0d0,183.0d0,
610
# 184.0d0,185.0d0,186.0d0,187.0d0,
611
# 188.0d0,189.0d0,190.0d0,195.0d0,
612
# 197.0d0,200.0d0,210.0d0,220.0d0,
613
# 230.0d0,240.0d0,250.0d0,260.0d0,
614
# 270.0d0,280.0d0,290.0d0,300.0d0,
615
# 310.0d0,320.0d0,330.0d0,335.0d0,
616
# 340.0d0,341.0d0,342.0d0,343.0d0,
617
# 344.0d0,345.0d0,345.3d0,345.5d0,
618
# 345.8d0,346.0d0,347.0d0,348.0d0,
619
# 349.0d0,350.0d0,351.0d0,352.0d0,
620
# 353.0d0,354.0d0,355.0d0,356.0d0,
621
# 357.0d0,358.0d0,359.0d0,360.0d0,
622
# 370.0d0,380.0d0,390.0d0,400.0d0,
623
# 410.0d0,420.0d0,430.0d0,440.0d0,
624
# 450.0d0,460.0d0,470.0d0,480.0d0,
625
# 490.0d0,500.0d0,510.0d0,520.0d0,
626
# 530.0d0,540.0d0,550.0d0,560.0d0,
627
# 570.0d0,580.0d0,590.0d0,600.0d0,
628
# 610.0d0,620.0d0,630.0d0,640.0d0,
629
# 650.0d0,660.0d0,670.0d0,680.0d0,
630
# 690.0d0,700.0d0,710.0d0,720.0d0,
631
# 730.0d0,740.0d0,750.0d0,760.0d0,
632
# 770.0d0,780.0d0,790.0d0,800.0d0,
633
# 810.0d0,820.0d0,830.0d0,840.0d0,
634
# 850.0d0,860.0d0,870.0d0,880.0d0,
635
# 890.0d0,900.0d0,910.0d0,920.0d0,
636
# 930.0d0,940.0d0,950.0d0,960.0d0,
637
# 970.0d0,980.0d0,990.0d0,1000.0d0/
638
*
639
DATA (yu(l),l=1,152)/4.187720723d0,4.549885500d0,
640
# 4.945324580d0,5.341899188d0,
641
# 5.713567160d0,5.874080782d0,5.983246138d0,5.991688858d0,
642
# 5.992041524d0,5.981547445d0,5.955788600d0,5.908389532d0,
643
# 5.829548098d0,5.703685183d0,5.508601342d0,5.222061322d0,
644
# 4.861934021d0,4.526068905d0,4.219144123d0,3.875133708d0,
645
# 3.522706643d0,3.196306140d0,2.906767625d0,2.652991136d0,
646
# 2.430066986d0,2.232861152d0,2.056219593d0,1.895596836d0,
647
# 1.747309568d0,1.607964744d0,1.474493488d0,1.342295851d0,
648
# 1.208919594d0,1.067309617d0,0.911812113d0,0.732879892d0,
649
# 0.522870581d0,0.290366470d0,0.068493581d0,-0.111577414d0,
650
# -0.279962552d0,-0.464332982d0,-0.654447961d0,-0.834336934d0,
651
# -0.999445210d0,-1.150871949d0,-1.287684579d0,-1.784688615d0,
652
# -1.924511090d0,-2.089081933d0,-2.402982024d0,-2.524611659d0,
653
# -2.525387652d0,-2.470417238d0,-2.399507117d0,-2.298953258d0,
654
# -2.192459462d0,-2.069700408d0,-1.973189437d0,-1.880904758d0,
655
# -1.807288585d0,-1.785086586d0,-1.826974607d0,-1.915990447d0,
656
# -2.062701208d0,-2.104793465d0,-2.152577512d0,-2.208540347d0,
657
# -2.276262386d0,-2.361739355d0,-2.391880716d0,-2.413272954d0,
658
# -2.447493257d0,-2.471820336d0,-2.622996766d0,-2.885017381d0,
659
# -3.696189261d0,-3.893536695d0,-4.011725974d0,-4.074946087d0,
660
# -4.137470131d0,-4.089349430d0,-4.136099918d0,-4.123934918d0,
661
# -4.136446777d0,-4.132255956d0,-4.118820717d0,-4.088009318d0,
662
# -3.788430916d0,-3.406565873d0,-3.050120268d0,-2.724817324d0,
663
# -2.410462214d0,-2.120474142d0,-1.845522666d0,-1.583137592d0,
664
# -1.337585290d0,-1.094144564d0,-0.897029424d0,-0.671551439d0,
665
# -0.478508120d0,-0.269184704d0,-0.072526771d0,0.115732816d0,
666
# 0.274371119d0,0.457987609d0,0.619147501d0,0.809852863d0,
667
# 0.968131131d0,1.121473977d0,1.290603896d0,1.462164252d0,
668
# 1.583836738d0,1.785815006d0,1.912473989d0,2.068564737d0,
669
# 2.235729083d0,2.398009567d0,2.538161792d0,2.681468958d0,
670
# 2.803535096d0,2.966400538d0,3.143956306d0,3.310124390d0,
671
# 3.429746138d0,3.580977910d0,3.727340256d0,3.890426351d0,
672
# 4.045983092d0,4.203320822d0,4.348209516d0,4.497115304d0,
673
# 4.654911988d0,4.823398533d0,4.983830013d0,5.127442851d0,
674
# 5.265056876d0,5.448819268d0,5.642707192d0,5.785441771d0,
675
# 5.936113103d0,6.108286154d0,6.258324908d0,6.415575321d0,
676
# 6.591607190d0,6.761124005d0,6.951304838d0,7.117437285d0,
677
# 7.283053264d0,7.453490877d0,7.637739897d0,7.799659928d0/
678
*
679
IF(top.eq.-1) THEN
680
n= 154
681
DO l=1,154
682
x(l)= xl(l)
683
y(l)= yl(l)
684
ENDDO
685
ELSEIF(top.eq.0) THEN
686
n= 151
687
DO l=1,151
688
x(l)= xc(l)
689
y(l)= yc(l)
690
ENDDO
691
ELSEIF(top.eq.1) THEN
692
n= 152
693
DO l=1,152
694
x(l)= xu(l)
695
y(l)= yu(l)
696
ENDDO
697
ENDIF
698
*
699
*.....First check if u is in the same interval found on the
700
* last call to Seval.............................................
701
*
702
IF ( (i.lt.1) .OR. (i.ge.n) ) i=1
703
IF ( (u.lt.x(i)) .OR. (u.ge.x(i+1)) ) THEN
704
i=1
705
*
706
* binary search
707
*
708
j= n+1
709
2468 k= (i+j)/2
710
IF (u.lt.x(k)) THEN
711
j= k
712
ELSE
713
i= k
714
ENDIF
715
IF (j.le.i+1) GOTO 2469
716
GOTO 2468
717
2469 CONTINUE
718
ENDIF
719
*
720
dx= u-x(i)
721
*
722
* evaluate the spline
723
*
724
f= y(i)+dx*(b(i)+dx*(c(i)+dx*d(i)))
725
fp= b(i)+dx*(TWO*c(i) + dx*THREE*d(i))
726
fpp= TWO*c(i) + dx*SIX*d(i)
727
fppp= SIX*d(i)
728
*
729
RETURN
730
*
731
END
732
733
*****************************************************************
734
C---electroweak corrections to H --> gamma gamma
735
*****************************************************************
736
737
double precision function gaga_elw(mtop,mh)
738
*
739
IMPLICIT NONE
740
*
741
INTEGER i,top,gdim
742
REAL*8 u,ui,value
743
REAL*8 vp,vpp,vppp
744
REAL*8 mtop,mh
745
REAL*8 x,x0,x1,x2,y0,y1,y2,a0,a1,a2
746
REAL*8 value0,value1,value2
747
REAL*8 bl(226),cl(226),dl(226)
748
REAL*8 bc(223),cc(223),dc(223)
749
REAL*8 bu(225),cu(225),du(225)
750
*
751
* u value of M_H at which the spline is to be evaluated
752
* top= -1,0,1 lower, central, upper value for m_top
753
*
754
u = mh
755
top = -1
756
gdim= 226
757
CALL FMMsplineSingle0(bl,cl,dl,top,gdim)
758
CALL Seval3Single0(u,bl,cl,dl,top,gdim,value0,vp,vpp,vppp)
759
top = 0
760
gdim= 223
761
CALL FMMsplineSingle0(bc,cc,dc,top,gdim)
762
CALL Seval3Single0(u,bc,cc,dc,top,gdim,value1,vp,vpp,vppp)
763
top = 1
764
gdim= 225
765
CALL FMMsplineSingle0(bu,cu,du,top,gdim)
766
CALL Seval3Single0(u,bu,cu,du,top,gdim,value2,vp,vpp,vppp)
767
x = mtop
768
x0 = 171.06d0
769
x1 = 172.64d0
770
x2 = 174.22d0
771
y0 = value0
772
y1 = value1
773
y2 = value2
774
A0=(X-X1)*(X-X2)/(X0-X1)/(X0-X2)
775
A1=(X-X0)*(X-X2)/(X1-X0)/(X1-X2)
776
A2=(X-X0)*(X-X1)/(X2-X0)/(X2-X1)
777
gaga_elw=(A0*Y0+A1*Y1+A2*Y2)/100.d0
778
c if(mh.gt.170.d0) gaga_elw = 0
779
RETURN
780
END
781
*
782
*-----------------------------------------------------------------------
783
*
784
c CONTAINS
785
*
786
SUBROUTINE FMMsplineSingle0(b,c,d,top,gdim)
787
*
788
* ---------------------------------------------------------------------------
789
*
790
* PURPOSE - Compute the coefficients b,c,d for a cubic interpolating spline
791
* so that the interpolated value is given by
792
* s(x) = y(k) + b(k)*(x-x(k)) + c(k)*(x-x(k))**2 + d(k)*(x-x(k))**3
793
* when x(k) <= x <= x(k+1)
794
* The end conditions match the third derivatives of the interpolated curve to
795
* the third derivatives of the unique polynomials thru the first four and
796
* last four points.
797
* Use Seval or Seval3 to evaluate the spline.
798
*
799
INTEGER k,n,i,top,gdim,l
800
*
801
REAL*8 xl(226),yl(226)
802
REAL*8 xc(223),yc(223)
803
REAL*8 xu(225),yu(225)
804
REAL*8 x(226),y(226)
805
*
806
REAL*8 b(gdim)
807
* linear coeff
808
*
809
REAL*8 c(gdim)
810
* quadratic coeff.
811
*
812
REAL*8 d(gdim)
813
* cubic coeff.
814
*
815
REAL*8 t
816
REAL*8 ZERO, TWO, THREE
817
PARAMETER(ZERO=0.0, TWO=2.0, THREE=3.0)
818
*
819
* The grid
820
*
821
DATA (xl(i),i=1,226)/100.0000d0,105.0000d0,110.0000d0,
822
# 115.0000d0,120.0000d0,125.0000d0,130.0000d0,135.0000d0,
823
# 140.0000d0,145.0000d0,150.0000d0,155.0000d0,156.0000d0,
824
# 157.0000d0,158.0000d0,159.0000d0,159.1000d0,159.2000d0,
825
# 159.3000d0,159.4000d0,159.5000d0,159.6000d0,159.7000d0,
826
# 159.8000d0,159.9000d0,160.0000d0,160.2500d0,160.5000d0,
827
# 160.7500d0,161.0000d0,161.5000d0,162.0000d0,162.5000d0,
828
# 163.0000d0,164.0000d0,165.0000d0,166.0000d0,167.0000d0,
829
# 168.0000d0,169.0000d0,170.0000d0,175.0000d0,179.0000d0,
830
# 180.0000d0,181.0000d0,182.0000d0,183.0000d0,184.0000d0,
831
# 185.0000d0,186.0000d0,187.0000d0,188.0000d0,189.0000d0,
832
# 190.0000d0,191.0000d0,192.0000d0,193.0000d0,194.0000d0,
833
# 195.0000d0,196.0000d0,197.0000d0,198.0000d0,199.0000d0,
834
# 200.0000d0,210.0000d0,220.0000d0,230.0000d0,240.0000d0,
835
# 250.0000d0,260.0000d0,270.0000d0,280.0000d0,290.0000d0,
836
# 300.0000d0,310.0000d0,320.0000d0,330.0000d0,340.0000d0,
837
# 340.5000d0,340.8000d0,341.0000d0,341.2000d0,341.4000d0,
838
# 341.6000d0,341.7000d0,341.8000d0,341.8500d0,341.9000d0,
839
# 341.9500d0,342.0000d0,342.0200d0,342.0400d0,342.0600d0,
840
# 342.0800d0,342.1000d0,342.1100d0,342.1120d0,342.1140d0,
841
# 342.1160d0,342.1180d0,342.1190d0,342.1196d0,342.1198d0,
842
# 342.1199d0,342.1204d0,342.1220d0,342.2000d0,342.3000d0,
843
# 342.4000d0,342.5000d0,342.6000d0,342.7000d0,342.8000d0,
844
# 342.9000d0,343.0000d0,343.2000d0,343.4000d0,343.6000d0,
845
# 343.8000d0,344.0000d0,345.0000d0,346.0000d0,347.0000d0,
846
# 348.0000d0,349.0000d0,350.0000d0,351.0000d0,352.0000d0,
847
# 353.0000d0,354.0000d0,355.0000d0,356.0000d0,357.0000d0,
848
# 358.0000d0,359.0000d0,360.0000d0,370.0000d0,380.0000d0,
849
# 390.0000d0,400.0000d0,410.0000d0,420.0000d0,430.0000d0,
850
# 440.0000d0,450.0000d0,460.0000d0,470.0000d0,480.0000d0,
851
# 490.0000d0,500.0000d0,510.0000d0,520.0000d0,530.0000d0,
852
# 540.0000d0,550.0000d0,560.0000d0,565.0000d0,570.0000d0,
853
# 575.0000d0,580.0000d0,585.0000d0,590.0000d0,595.0000d0,
854
# 600.0000d0,605.0000d0,610.0000d0,615.0000d0,620.0000d0,
855
# 625.0000d0,630.0000d0,635.0000d0,640.0000d0,641.0000d0,
856
# 642.0000d0,643.0000d0,644.0000d0,645.0000d0,646.0000d0,
857
# 647.0000d0,648.0000d0,649.0000d0,650.0000d0,651.0000d0,
858
# 652.0000d0,653.0000d0,654.0000d0,655.0000d0,656.0000d0,
859
# 657.0000d0,658.0000d0,659.0000d0,660.0000d0,670.0000d0,
860
# 680.0000d0,690.0000d0,700.0000d0,710.0000d0,720.0000d0,
861
# 730.0000d0,740.0000d0,750.0000d0,760.0000d0,770.0000d0,
862
# 780.0000d0,790.0000d0,800.0000d0,810.0000d0,820.0000d0,
863
# 830.0000d0,840.0000d0,850.0000d0,860.0000d0,870.0000d0,
864
# 880.0000d0,890.0000d0,900.0000d0,910.0000d0,920.0000d0,
865
# 930.0000d0,940.0000d0,950.0000d0,960.0000d0,970.0000d0,
866
# 980.0000d0,990.0000d0,1000.0000d0/
867
*
868
DATA (yl(i),i=1,226)/-2.927515523d0,-2.713875277d0,-2.477263056d0,
869
# -2.215702562d0, -1.926872830d0, -1.607794154d0, -1.254642316d0,
870
# -0.862759386d0, -0.424192821d0, 0.071393692d0, 0.642180822d0,
871
# 1.311566449d0, 1.458302972d0, 1.610089467d0, 1.771001228d0,
872
# 1.958842368d0, 1.980856139d0, 2.003769992d0, 2.027683661d0,
873
# 2.052701483d0, 2.078930176d0, 2.106476183d0, 2.135441928d0,
874
# 2.165921345d0, 2.197994437d0, 2.231721476d0, 2.323417121d0,
875
# 2.425109533d0, 2.534671123d0, 2.648366651d0, 2.869451789d0,
876
# 3.056555683d0, 3.197055532d0, 3.294656543d0, 3.399442397d0,
877
# 3.438265245d0, 3.450036365d0, 3.449871888d0, 3.442815310d0,
878
# 3.428262691d0, 3.411797242d0, 3.363961271d0, 3.295179282d0,
879
# 3.287861332d0, 3.325219720d0, 3.470009593d0, 3.737173473d0,
880
# 4.021846379d0, 4.241836705d0, 4.394931599d0, 4.501421364d0,
881
# 4.577551451d0, 4.633363121d0, 4.675411596d0, 4.707771996d0,
882
# 4.732593619d0, 4.751501682d0, 4.766245197d0, 4.778023555d0,
883
# 4.787199047d0, 4.793849566d0, 4.798265576d0, 4.801110206d0,
884
# 4.802711433d0, 4.767182955d0, 4.677220707d0, 4.558671105d0,
885
# 4.423406368d0, 4.276574444d0, 4.119375511d0, 3.974083777d0,
886
# 3.819733346d0, 3.671217773d0, 3.517370222d0, 3.344189966d0,
887
# 3.116802834d0, 2.754066836d0, 1.809839704d0, 1.691727991d0,
888
# 1.607871864d0, 1.544493194d0, 1.473234377d0, 1.390992309d0,
889
# 1.292466985d0, 1.234537350d0, 1.167169430d0, 1.128863846d0,
890
# 1.085963030d0, 1.037416400d0, 0.979461824d0, 0.952439442d0,
891
# 0.922628624d0, 0.888386832d0, 0.846809865d0, 0.791742524d0,
892
# 0.752347274d0, 0.742175952d0, 0.730422928d0, 0.716838887d0,
893
# 0.697956938d0, 0.684261787d0, 0.672163545d0, 0.665924128d0,
894
# 0.661453423d0, 0.652439216d0, 0.656998002d0, 0.677275168d0,
895
# 0.699459467d0, 0.721833458d0, 0.738243489d0, 0.751895487d0,
896
# 0.769698326d0, 0.779779847d0, 0.799512664d0, 0.816893187d0,
897
# 0.842794338d0, 0.874519045d0, 0.916496501d0, 0.950264630d0,
898
# 0.968728477d0, 1.105576625d0, 1.278360270d0, 1.391736231d0,
899
# 1.542276157d0, 1.675869257d0, 1.795595314d0, 1.923303286d0,
900
# 2.025352084d0, 2.157383377d0, 2.260417650d0, 2.361640519d0,
901
# 2.479888170d0, 2.572642417d0, 2.688809606d0, 2.796999182d0,
902
# 2.879392617d0, 3.848596335d0, 4.732837134d0, 5.525766484d0,
903
# 6.243989288d0, 6.975830059d0, 7.607337162d0, 8.167782071d0,
904
# 8.790374186d0, 9.422028844d0, 9.922844338d0, 10.426856521d0,
905
# 10.933372206d0, 11.329434560d0, 11.735173309d0, 11.795158880d0,
906
# 11.720854043d0, 11.455441492d0, 10.516061928d0, 9.017310079d0,
907
# 7.000279246d0, 5.560369962d0, 3.715234200d0, 1.373399680d0,
908
# -1.698051464d0, -5.394210410d0, -9.477907268d0,-13.677692505d0,
909
# -18.256012300d0,-24.029072069d0,-29.992304980d0,-35.639341296d0,
910
# -41.517060283d0,-46.430415213d0,-50.471411203d0,-53.247623713d0,
911
# -55.184140125d0,-55.523807703d0,-55.749322940d0,-55.785869793d0,
912
# -56.086101947d0,-56.289055861d0,-56.369321707d0,-56.555818993d0,
913
# -56.632365789d0,-56.616832544d0,-56.709629986d0,-56.539655434d0,
914
# -56.366567734d0,-56.272247546d0,-56.170833793d0,-55.984866548d0,
915
# -55.776968735d0,-55.354292076d0,-55.136594171d0,-54.856130174d0,
916
# -54.390381982d0,-50.585126291d0,-46.161103219d0,-41.187366158d0,
917
# -36.554683053d0,-32.185245783d0,-28.097939638d0,-24.608441494d0,
918
# -21.643195994d0,-18.984106588d0,-16.998672281d0,-15.273879570d0,
919
# -13.373334219d0,-11.294832079d0, -9.584938331d0, -8.286389135d0,
920
# -6.838755740d0, -5.409090746d0, -4.588339018d0, -3.705250749d0,
921
# -2.476432078d0, -1.367536281d0, -0.414988031d0, 0.510347996d0,
922
# 1.376925204d0, 2.276754999d0, 3.203065379d0, 3.629617497d0,
923
# 3.821185286d0, 4.122178708d0, 4.804270233d0, 5.929961061d0,
924
# 7.072311664d0, 7.458878393d0, 7.823110180d0/
925
*
926
DATA (xc(i),i=1,223)/100.0000d0,105.0000d0,110.0000d0,
927
# 115.0000d0,120.0000d0,125.0000d0,130.0000d0,135.0000d0,
928
# 140.0000d0,145.0000d0,150.0000d0,155.0000d0,156.0000d0,
929
# 157.0000d0,158.0000d0,159.0000d0,159.1000d0,159.2000d0,
930
# 159.3000d0,159.4000d0,159.5000d0,159.6000d0,159.7000d0,
931
# 159.8000d0,159.9000d0,160.0000d0,160.2500d0,160.5000d0,
932
# 160.7500d0,161.0000d0,161.5000d0,162.0000d0,162.5000d0,
933
# 163.0000d0,164.0000d0,165.0000d0,166.0000d0,167.0000d0,
934
# 168.0000d0,169.0000d0,170.0000d0,175.0000d0,178.5000d0,
935
# 179.0000d0,179.5000d0,180.0000d0,180.5000d0,181.0000d0,
936
# 181.5000d0,182.0000d0,182.5000d0,183.0000d0,184.0000d0,
937
# 185.0000d0,186.0000d0,187.0000d0,188.0000d0,189.0000d0,
938
# 190.0000d0,191.0000d0,192.0000d0,193.0000d0,194.0000d0,
939
# 195.0000d0,196.0000d0,197.0000d0,198.0000d0,199.0000d0,
940
# 200.0000d0,210.0000d0,220.0000d0,230.0000d0,240.0000d0,
941
# 250.0000d0,260.0000d0,270.0000d0,280.0000d0,290.0000d0,
942
# 300.0000d0,310.0000d0,320.0000d0,330.0000d0,340.0000d0,
943
# 341.0000d0,342.0000d0,343.0000d0,344.0000d0,344.2000d0,
944
# 344.4000d0,344.6000d0,344.8000d0,344.9000d0,345.0000d0,
945
# 345.0500d0,345.1000d0,345.1500d0,345.2000d0,345.2200d0,
946
# 345.2400d0,345.2600d0,345.2700d0,345.2720d0,345.2740d0,
947
# 345.2760d0,345.2780d0,345.2790d0,345.2796d0,345.2798d0,
948
# 345.2799d0,345.2804d0,345.2820d0,345.2900d0,345.3500d0,
949
# 345.4500d0,345.5200d0,345.6000d0,345.7000d0,345.8000d0,
950
# 346.0000d0,347.0000d0,348.0000d0,349.0000d0,350.0000d0,
951
# 351.0000d0,352.0000d0,353.0000d0,354.0000d0,355.0000d0,
952
# 356.0000d0,357.0000d0,358.0000d0,359.0000d0,360.0000d0,
953
# 370.0000d0,380.0000d0,390.0000d0,400.0000d0,410.0000d0,
954
# 420.0000d0,430.0000d0,440.0000d0,450.0000d0,460.0000d0,
955
# 470.0000d0,480.0000d0,490.0000d0,500.0000d0,510.0000d0,
956
# 520.0000d0,530.0000d0,540.0000d0,550.0000d0,560.0000d0,
957
# 565.0000d0,570.0000d0,575.0000d0,580.0000d0,585.0000d0,
958
# 590.0000d0,595.0000d0,600.0000d0,605.0000d0,610.0000d0,
959
# 615.0000d0,620.0000d0,625.0000d0,630.0000d0,635.0000d0,
960
# 640.0000d0,641.0000d0,642.0000d0,643.0000d0,644.0000d0,
961
# 645.0000d0,646.0000d0,647.0000d0,648.0000d0,649.0000d0,
962
# 650.0000d0,651.0000d0,652.0000d0,653.0000d0,654.0000d0,
963
# 655.0000d0,656.0000d0,657.0000d0,658.0000d0,659.0000d0,
964
# 660.0000d0,670.0000d0,680.0000d0,690.0000d0,700.0000d0,
965
# 710.0000d0,720.0000d0,730.0000d0,740.0000d0,750.0000d0,
966
# 760.0000d0,770.0000d0,780.0000d0,790.0000d0,800.0000d0,
967
# 810.0000d0,820.0000d0,830.0000d0,840.0000d0,850.0000d0,
968
# 860.0000d0,870.0000d0,880.0000d0,890.0000d0,900.0000d0,
969
# 910.0000d0,920.0000d0,930.0000d0,940.0000d0,950.0000d0,
970
# 960.0000d0,970.0000d0,980.0000d0,990.0000d0,1000.0000d0/
971
*
972
DATA (yc(i),i=1,223)/-2.969458664d0,-2.755568498d0,-2.518689112d0,
973
# -2.256840599d0, -1.967703151d0, -1.648299088d0, -1.294794285d0,
974
# -0.902525232d0, -0.463541061d0, 0.032517848d0, 0.603841958d0,
975
# 1.273887882d0, 1.420789352d0, 1.572747980d0, 1.733843353d0,
976
# 1.921879203d0, 1.943912705d0, 1.966846336d0, 1.990779880d0,
977
# 2.015817728d0, 2.042066621d0, 2.069633010d0, 2.098619322d0,
978
# 2.129119489d0, 2.161213532d0, 2.194961744d0, 2.286711412d0,
979
# 2.388459445d0, 2.498078035d0, 2.611830784d0, 2.833024715d0,
980
# 3.020227220d0, 3.160820150d0, 3.258502948d0, 3.363381747d0,
981
# 3.402305539d0, 3.414208425d0, 3.414049050d0, 3.407053730d0,
982
# 3.392629803d0, 3.376195597d0, 3.328499111d0, 3.268651902d0,
983
# 3.259860594d0, 3.253486002d0, 3.252546948d0, 3.262154622d0,
984
# 3.289924633d0, 3.345085115d0, 3.434732330d0, 3.557729513d0,
985
# 3.701897111d0, 3.986569984d0, 4.206557658d0, 4.359682402d0,
986
# 4.466176116d0, 4.542300703d0, 4.598184529d0, 4.640273713d0,
987
# 4.672604523d0, 4.697438255d0, 4.716361927d0, 4.731103154d0,
988
# 4.742925539d0, 4.752092918d0, 4.758670785d0, 4.763130070d0,
989
# 4.766074843d0, 4.767748196d0, 4.732269356d0, 4.642413616d0,
990
# 4.524181528d0, 4.389510854d0, 4.242822112d0, 4.086217524d0,
991
# 3.941680893d0, 3.789691168d0, 3.644241595d0, 3.494351920d0,
992
# 3.331867639d0, 3.122929519d0, 2.816374576d0, 2.186038507d0,
993
# 2.068538445d0, 1.926511260d0, 1.746359564d0, 1.494945494d0,
994
# 1.428923128d0, 1.354291504d0, 1.267892764d0, 1.163289817d0,
995
# 1.100587592d0, 1.027248191d0, 0.984553881d0, 0.936444199d0,
996
# 0.879742362d0, 0.808975906d0, 0.774145495d0, 0.731724792d0,
997
# 0.675557165d0, 0.635406715d0, 0.625004410d0, 0.613030568d0,
998
# 0.599164524d0, 0.579800865d0, 0.565857087d0, 0.553503488d0,
999
# 0.547057831d0, 0.542411799d0, 0.533743538d0, 0.538671942d0,
1000
# 0.542814415d0, 0.562187598d0, 0.575771207d0, 0.586370628d0,
1001
# 0.605546923d0, 0.622661850d0, 0.639021282d0, 0.683852301d0,
1002
# 0.846885709d0, 0.995547199d0, 1.139428677d0, 1.270752033d0,
1003
# 1.403862461d0, 1.540060220d0, 1.676615197d0, 1.791154110d0,
1004
# 1.909977528d0, 2.023870438d0, 2.138393119d0, 2.233199625d0,
1005
# 2.358706648d0, 2.446271802d0, 3.459751114d0, 4.411864477d0,
1006
# 5.220622691d0, 5.988967245d0, 6.706065217d0, 7.359912016d0,
1007
# 7.993290467d0, 8.675262256d0, 9.327550412d0, 9.917129519d0,
1008
# 10.373037136d0, 11.009210530d0, 11.516244354d0, 11.985967329d0,
1009
# 12.198319618d0, 12.357011718d0, 12.311558072d0, 11.754465713d0,
1010
# 10.617092695d0, 9.043395280d0, 8.002055503d0, 6.601203443d0,
1011
# 4.587232698d0, 1.906399387d0, -1.281177539d0, -5.052281649d0,
1012
# -9.213421228d0,-13.868726558d0,-19.489194521d0,-25.988103189d0,
1013
# -32.734525726d0,-39.402198052d0,-45.903352541d0,-51.189352305d0,
1014
# -55.322071720d0,-58.186905136d0,-58.653812374d0,-59.048160142d0,
1015
# -59.279596622d0,-59.992120796d0,-60.271154555d0,-60.549383981d0,
1016
# -60.729563975d0,-60.954540486d0,-61.084961955d0,-61.164811836d0,
1017
# -61.186911575d0,-61.109847454d0,-61.108201824d0,-60.992397754d0,
1018
# -60.768938455d0,-60.645603901d0,-60.224362756d0,-59.860911894d0,
1019
# -59.548720101d0,-59.301307840d0,-55.034929471d0,-49.745997496d0,
1020
# -44.181440809d0,-39.061582838d0,-34.158470096d0,-29.603188102d0,
1021
# -25.799341483d0,-22.634533844d0,-19.804942358d0,-17.654326276d0,
1022
# -15.913630982d0,-13.830699918d0,-11.633319645d0,-10.006868842d0,
1023
# -8.532937676d0, -7.132846889d0, -5.627566626d0, -4.758444504d0,
1024
# -3.887401469d0, -2.729580751d0, -1.512994295d0, -0.550651909d0,
1025
# 0.355022046d0, 1.248689020d0, 2.123963438d0, 3.014801413d0,
1026
# 3.484406765d0, 3.696913629d0, 4.008426964d0, 4.651255618d0,
1027
# 5.840647264d0, 6.935940541d0, 7.364646872d0, 7.670905366d0/
1028
*
1029
DATA (xu(i),i=1,225)/100.0000d0,105.0000d0,110.0000d0,
1030
# 115.0000d0,120.0000d0,125.0000d0,130.0000d0,135.0000d0,
1031
# 140.0000d0,145.0000d0,150.0000d0,155.0000d0,156.0000d0,
1032
# 157.0000d0,158.0000d0,159.0000d0,159.1000d0,159.2000d0,
1033
# 159.3000d0,159.4000d0,159.5000d0,159.6000d0,159.7000d0,
1034
# 159.8000d0,159.9000d0,160.0000d0,160.2500d0,160.5000d0,
1035
# 160.7500d0,161.0000d0,161.5000d0,162.0000d0,162.5000d0,
1036
# 163.0000d0,164.0000d0,165.0000d0,166.0000d0,167.0000d0,
1037
# 168.0000d0,169.0000d0,170.0000d0,175.0000d0,178.5000d0,
1038
# 179.0000d0,179.5000d0,180.0000d0,180.5000d0,181.0000d0,
1039
# 181.5000d0,182.0000d0,182.5000d0,183.0000d0,184.0000d0,
1040
# 185.0000d0,186.0000d0,187.0000d0,188.0000d0,189.0000d0,
1041
# 190.0000d0,191.0000d0,192.0000d0,193.0000d0,194.0000d0,
1042
# 195.0000d0,196.0000d0,197.0000d0,198.0000d0,199.0000d0,
1043
# 200.0000d0,210.0000d0,220.0000d0,230.0000d0,240.0000d0,
1044
# 250.0000d0,260.0000d0,270.0000d0,280.0000d0,290.0000d0,
1045
# 300.0000d0,310.0000d0,320.0000d0,330.0000d0,340.0000d0,
1046
# 341.0000d0,342.0000d0,343.0000d0,344.0000d0,345.0000d0,
1047
# 345.5000d0,345.8000d0,346.0000d0,347.0000d0,347.2000d0,
1048
# 347.4000d0,347.6000d0,347.8000d0,348.0000d0,348.1000d0,
1049
# 348.2000d0,348.2500d0,348.3000d0,348.3200d0,348.3400d0,
1050
# 348.3600d0,348.3800d0,348.4000d0,348.4100d0,348.4200d0,
1051
# 348.4300d0,348.4320d0,348.4340d0,348.4360d0,348.4380d0,
1052
# 348.4390d0,348.4396d0,348.4398d0,348.4399d0,348.4404d0,
1053
# 348.4420d0,348.4500d0,348.5000d0,348.7000d0,349.0000d0,
1054
# 349.4000d0,350.0000d0,351.0000d0,352.0000d0,353.0000d0,
1055
# 354.0000d0,355.0000d0,356.0000d0,357.0000d0,358.0000d0,
1056
# 359.0000d0,360.0000d0,370.0000d0,380.0000d0,390.0000d0,
1057
# 400.0000d0,410.0000d0,420.0000d0,430.0000d0,440.0000d0,
1058
# 450.0000d0,460.0000d0,470.0000d0,480.0000d0,490.0000d0,
1059
# 500.0000d0,510.0000d0,520.0000d0,530.0000d0,540.0000d0,
1060
# 550.0000d0,560.0000d0,565.0000d0,570.0000d0,575.0000d0,
1061
# 580.0000d0,585.0000d0,590.0000d0,595.0000d0,600.0000d0,
1062
# 605.0000d0,610.0000d0,615.0000d0,620.0000d0,625.0000d0,
1063
# 630.0000d0,635.0000d0,640.0000d0,641.0000d0,642.0000d0,
1064
# 643.0000d0,644.0000d0,645.0000d0,646.0000d0,647.0000d0,
1065
# 648.0000d0,649.0000d0,650.0000d0,651.0000d0,652.0000d0,
1066
# 653.0000d0,654.0000d0,655.0000d0,656.0000d0,657.0000d0,
1067
# 658.0000d0,659.0000d0,660.0000d0,670.0000d0,680.0000d0,
1068
# 690.0000d0,700.0000d0,710.0000d0,720.0000d0,730.0000d0,
1069
# 740.0000d0,750.0000d0,760.0000d0,770.0000d0,780.0000d0,
1070
# 790.0000d0,800.0000d0,810.0000d0,820.0000d0,830.0000d0,
1071
# 840.0000d0,850.0000d0,860.0000d0,870.0000d0,880.0000d0,
1072
# 890.0000d0,900.0000d0,910.0000d0,920.0000d0,930.0000d0,
1073
# 940.0000d0,950.0000d0,960.0000d0,970.0000d0,980.0000d0,
1074
# 990.0000d0,1000.0000d0/
1075
*
1076
DATA (yu(i),i=1,225)/-3.011799533d0,-2.797658814d0,-2.560511342d0,
1077
# -2.298374579d0, -2.008927458d0, -1.689194746d0, -1.335334782d0,
1078
# -0.942678003d0, -0.503267763d0, -0.006741983d0, 0.565124142d0,
1079
# 1.235840116d0, 1.382905273d0, 1.535044791d0, 1.696328604d0,
1080
# 1.884567714d0, 1.906622347d0, 1.929577239d0, 1.953532153d0,
1081
# 1.978591436d0, 2.004861792d0, 2.032449634d0, 2.061457355d0,
1082
# 2.091978868d0, 2.124094156d0, 2.157863499d0, 2.249665540d0,
1083
# 2.351465763d0, 2.461136573d0, 2.574941494d0, 2.796240093d0,
1084
# 2.983545301d0, 3.124228166d0, 3.221982342d0, 3.326972382d0,
1085
# 3.365993805d0, 3.378023193d0, 3.377887062d0, 3.370931083d0,
1086
# 3.356634680d0, 3.340265619d0, 3.292682936d0, 3.232961811d0,
1087
# 3.224194020d0, 3.217828260d0, 3.216889093d0, 3.226496714d0,
1088
# 3.254271720d0, 3.309440914d0, 3.399097855d0, 3.522103315d0,
1089
# 3.666276367d0, 3.950949463d0, 4.170926862d0, 4.324070476d0,
1090
# 4.430588837d0, 4.506706461d0, 4.562622634d0, 4.604788444d0,
1091
# 4.637112788d0, 4.661928180d0, 4.680880990d0, 4.695625131d0,
1092
# 4.707459275d0, 4.716670603d0, 4.723178694d0, 4.727608289d0,
1093
# 4.730650318d0, 4.732442748d0, 4.697030978d0, 4.607168713d0,
1094
# 4.489397175d0, 4.355888392d0, 4.208987612d0, 4.052707904d0,
1095
# 3.908019896d0, 3.759968739d0, 3.615348691d0, 3.469932507d0,
1096
# 3.315332978d0, 3.122286809d0, 2.855004727d0, 2.377722123d0,
1097
# 2.300655190d0, 2.213908305d0, 2.115631830d0, 2.000162289d0,
1098
# 1.860840425d0, 1.779119559d0, 1.724912832d0, 1.685992661d0,
1099
# 1.446048759d0, 1.384730617d0, 1.316264190d0, 1.238304389d0,
1100
# 1.147571101d0, 1.036738128d0, 0.969143795d0, 0.888640045d0,
1101
# 0.841031364d0, 0.785526212d0, 0.759880178d0, 0.731970096d0,
1102
# 0.701136254d0, 0.665654946d0, 0.622492994d0, 0.596679361d0,
1103
# 0.565315972d0, 0.524528227d0, 0.513928834d0, 0.501826683d0,
1104
# 0.487735928d0, 0.468019598d0, 0.454013349d0, 0.441470198d0,
1105
# 0.434889761d0, 0.430148561d0, 0.421464835d0, 0.429200773d0,
1106
# 0.429255984d0, 0.437861500d0, 0.480919141d0, 0.544499260d0,
1107
# 0.606035270d0, 0.710025546d0, 0.864539937d0, 1.009890014d0,
1108
# 1.156589498d0, 1.268264088d0, 1.408793136d0, 1.528293889d0,
1109
# 1.642389803d0, 1.771929787d0, 1.875426734d0, 1.979445634d0,
1110
# 3.092059170d0, 4.047318673d0, 4.937427406d0, 5.684883866d0,
1111
# 6.473364254d0, 7.141842873d0, 7.823212681d0, 8.484967861d0,
1112
# 9.214907909d0, 9.784981877d0, 10.383431544d0, 11.054527450d0,
1113
# 11.588357911d0, 12.165403348d0, 12.507103232d0, 12.848157153d0,
1114
# 12.975158246d0, 12.770167330d0, 11.991429721d0, 11.123821384d0,
1115
# 10.241165300d0, 9.157568585d0, 7.611314332d0, 5.499567320d0,
1116
# 2.603705830d0, -0.579201122d0, -4.348429690d0, -8.823492981d0,
1117
# -14.416229606d0,-20.848895478d0,-28.011934618d0,-36.077657061d0,
1118
# -43.726529167d0,-50.682419374d0,-56.593917330d0,-60.848713817d0,
1119
# -61.464423821d0,-62.283661600d0,-62.930084228d0,-63.549395305d0,
1120
# -64.133867221d0,-64.465495513d0,-64.909336911d0,-65.479678788d0,
1121
# -65.703953492d0,-66.113782243d0,-66.222507344d0,-66.278786408d0,
1122
# -66.393237110d0,-66.366867383d0,-66.202199439d0,-66.029041364d0,
1123
# -65.774112798d0,-65.485847916d0,-65.247943137d0,-64.801055924d0,
1124
# -60.101647018d0,-54.084245796d0,-47.571709079d0,-41.661629810d0,
1125
# -36.168693788d0,-31.212363345d0,-27.000704767d0,-23.605327875d0,
1126
# -20.580041179d0,-18.274370612d0,-16.434309906d0,-14.327218251d0,
1127
# -12.021979554d0,-10.331796926d0, -8.883031936d0, -7.348435651d0,
1128
# -5.780362152d0, -4.953756827d0, -4.070387325d0, -2.840184650d0,
1129
# -1.678475825d0, -0.679391004d0, 0.223829200d0, 1.086204496d0,
1130
# 2.026633865d0, 2.923942775d0, 3.306978920d0, 3.579966278d0,
1131
# 3.863101723d0, 4.490999475d0, 5.727923150d0, 6.829991246d0,
1132
# 7.212458092d0, 7.551889308d0/
1133
*
1134
IF(top.eq.-1) THEN
1135
n= 226
1136
DO l=1,226
1137
x(l)= xl(l)
1138
y(l)= yl(l)
1139
c write(35,*)l,x(l),y(l)
1140
ENDDO
1141
ELSEIF(top.eq.0) THEN
1142
n= 223
1143
DO l=1,223
1144
x(l)= xc(l)
1145
y(l)= yc(l)
1146
c write(36,*)l,x(l),y(l)
1147
ENDDO
1148
ELSEIF(top.eq.1) THEN
1149
n= 225
1150
DO l=1,225
1151
x(l)= xu(l)
1152
y(l)= yu(l)
1153
c write(37,*)l,x(l),y(l)
1154
ENDDO
1155
ENDIF
1156
1157
*.....Set up tridiagonal system.........................................
1158
* b=diagonal, d=offdiagonal, c=right-hand side
1159
*
1160
d(1)= x(2)-x(1)
1161
c(2)= (y(2)-y(1))/d(1)
1162
DO k= 2,n-1
1163
d(k)= x(k+1)-x(k)
1164
b(k)= TWO*(d(k-1)+d(k))
1165
c(k+1)= (y(k+1)-y(k))/d(k)
1166
c(k)= c(k+1)-c(k)
1167
END DO
1168
*
1169
*.....End conditions. third derivatives at x(1) and x(n) obtained
1170
* from divided differences.......................................
1171
*
1172
b(1)= -d(1)
1173
b(n)= -d(n-1)
1174
c(1)= ZERO
1175
c(n)= ZERO
1176
IF (n.gt.3) THEN
1177
c(1)= c(3)/(x(4)-x(2))-c(2)/(x(3)-x(1))
1178
c(n)= c(n-1)/(x(n)-x(n-2))-c(n-2)/(x(n-1)-x(n-3))
1179
c(1)= c(1)*d(1)*d(1)/(x(4)-x(1))
1180
c(n)= -c(n)*d(n-1)*d(n-1)/(x(n)-x(n-3))
1181
END IF
1182
*
1183
DO k=2,n ! forward elimination
1184
t= d(k-1)/b(k-1)
1185
b(k)= b(k)-t*d(k-1)
1186
c(k)= c(k)-t*c(k-1)
1187
END DO
1188
*
1189
c(n)= c(n)/b(n)
1190
*
1191
* back substitution ( makes c the sigma of text)
1192
*
1193
DO k=n-1,1,-1
1194
c(k)= (c(k)-d(k)*c(k+1))/b(k)
1195
END DO
1196
*
1197
*.....Compute polynomial coefficients...................................
1198
*
1199
b(n)= (y(n)-y(n-1))/d(n-1)+d(n-1)*(c(n-1)+c(n)+c(n))
1200
DO k=1,n-1
1201
b(k)= (y(k+1)-y(k))/d(k)-d(k)*(c(k+1)+c(k)+c(k))
1202
d(k)= (c(k+1)-c(k))/d(k)
1203
c(k)= THREE*c(k)
1204
END DO
1205
c(n)= THREE*c(n)
1206
d(n)= d(n-1)
1207
*
1208
RETURN
1209
*
1210
END
1211
*
1212
*------------------------------------------------------------------------
1213
*
1214
SUBROUTINE Seval3Single0(u,b,c,d,top,gdim,f,fp,fpp,fppp)
1215
*
1216
* ---------------------------------------------------------------------------
1217
*
1218
* PURPOSE - Evaluate the cubic spline function
1219
* Seval=y(i)+b(i)*(u-x(i))+c(i)*(u-x(i))**2+d(i)*(u-x(i))**3
1220
* where x(i) <= u < x(i+1)
1221
*
1222
* NOTES- if u<x(1), i=1 is used;if u>x(n), i=n is used
1223
*
1224
REAL*8 u
1225
* abscissa at which the spline is to be evaluated
1226
*
1227
INTEGER j,k,n,l,top,gdim
1228
*
1229
REAL*8 xl(226),yl(226)
1230
REAL*8 xc(223),yc(223)
1231
REAL*8 xu(225),yu(225)
1232
REAL*8 x(226),y(226)
1233
REAL*8 b(gdim),c(gdim),d(gdim)
1234
* linear,quadratic,cubic coeff
1235
*
1236
REAL*8 f,fp,fpp,fppp
1237
* function, 1st,2nd,3rd deriv
1238
*
1239
INTEGER i
1240
REAL*8 dx
1241
REAL*8 ZERO, TWO, THREE
1242
PARAMETER(ZERO=0.0, TWO=2.0, THREE=3.0)
1243
SAVE i
1244
DATA i/1/
1245
*
1246
* The grid
1247
*
1248
DATA (xl(i),i=1,226)/100.0000d0,105.0000d0,110.0000d0,
1249
# 115.0000d0,120.0000d0,125.0000d0,130.0000d0,135.0000d0,
1250
# 140.0000d0,145.0000d0,150.0000d0,155.0000d0,156.0000d0,
1251
# 157.0000d0,158.0000d0,159.0000d0,159.1000d0,159.2000d0,
1252
# 159.3000d0,159.4000d0,159.5000d0,159.6000d0,159.7000d0,
1253
# 159.8000d0,159.9000d0,160.0000d0,160.2500d0,160.5000d0,
1254
# 160.7500d0,161.0000d0,161.5000d0,162.0000d0,162.5000d0,
1255
# 163.0000d0,164.0000d0,165.0000d0,166.0000d0,167.0000d0,
1256
# 168.0000d0,169.0000d0,170.0000d0,175.0000d0,179.0000d0,
1257
# 180.0000d0,181.0000d0,182.0000d0,183.0000d0,184.0000d0,
1258
# 185.0000d0,186.0000d0,187.0000d0,188.0000d0,189.0000d0,
1259
# 190.0000d0,191.0000d0,192.0000d0,193.0000d0,194.0000d0,
1260
# 195.0000d0,196.0000d0,197.0000d0,198.0000d0,199.0000d0,
1261
# 200.0000d0,210.0000d0,220.0000d0,230.0000d0,240.0000d0,
1262
# 250.0000d0,260.0000d0,270.0000d0,280.0000d0,290.0000d0,
1263
# 300.0000d0,310.0000d0,320.0000d0,330.0000d0,340.0000d0,
1264
# 340.5000d0,340.8000d0,341.0000d0,341.2000d0,341.4000d0,
1265
# 341.6000d0,341.7000d0,341.8000d0,341.8500d0,341.9000d0,
1266
# 341.9500d0,342.0000d0,342.0200d0,342.0400d0,342.0600d0,
1267
# 342.0800d0,342.1000d0,342.1100d0,342.1120d0,342.1140d0,
1268
# 342.1160d0,342.1180d0,342.1190d0,342.1196d0,342.1198d0,
1269
# 342.1199d0,342.1204d0,342.1220d0,342.2000d0,342.3000d0,
1270
# 342.4000d0,342.5000d0,342.6000d0,342.7000d0,342.8000d0,
1271
# 342.9000d0,343.0000d0,343.2000d0,343.4000d0,343.6000d0,
1272
# 343.8000d0,344.0000d0,345.0000d0,346.0000d0,347.0000d0,
1273
# 348.0000d0,349.0000d0,350.0000d0,351.0000d0,352.0000d0,
1274
# 353.0000d0,354.0000d0,355.0000d0,356.0000d0,357.0000d0,
1275
# 358.0000d0,359.0000d0,360.0000d0,370.0000d0,380.0000d0,
1276
# 390.0000d0,400.0000d0,410.0000d0,420.0000d0,430.0000d0,
1277
# 440.0000d0,450.0000d0,460.0000d0,470.0000d0,480.0000d0,
1278
# 490.0000d0,500.0000d0,510.0000d0,520.0000d0,530.0000d0,
1279
# 540.0000d0,550.0000d0,560.0000d0,565.0000d0,570.0000d0,
1280
# 575.0000d0,580.0000d0,585.0000d0,590.0000d0,595.0000d0,
1281
# 600.0000d0,605.0000d0,610.0000d0,615.0000d0,620.0000d0,
1282
# 625.0000d0,630.0000d0,635.0000d0,640.0000d0,641.0000d0,
1283
# 642.0000d0,643.0000d0,644.0000d0,645.0000d0,646.0000d0,
1284
# 647.0000d0,648.0000d0,649.0000d0,650.0000d0,651.0000d0,
1285
# 652.0000d0,653.0000d0,654.0000d0,655.0000d0,656.0000d0,
1286
# 657.0000d0,658.0000d0,659.0000d0,660.0000d0,670.0000d0,
1287
# 680.0000d0,690.0000d0,700.0000d0,710.0000d0,720.0000d0,
1288
# 730.0000d0,740.0000d0,750.0000d0,760.0000d0,770.0000d0,
1289
# 780.0000d0,790.0000d0,800.0000d0,810.0000d0,820.0000d0,
1290
# 830.0000d0,840.0000d0,850.0000d0,860.0000d0,870.0000d0,
1291
# 880.0000d0,890.0000d0,900.0000d0,910.0000d0,920.0000d0,
1292
# 930.0000d0,940.0000d0,950.0000d0,960.0000d0,970.0000d0,
1293
# 980.0000d0,990.0000d0,1000.0000d0/
1294
*
1295
DATA (yl(i),i=1,226)/-2.927515523d0,-2.713875277d0,-2.477263056d0,
1296
# -2.215702562d0, -1.926872830d0, -1.607794154d0, -1.254642316d0,
1297
# -0.862759386d0, -0.424192821d0, 0.071393692d0, 0.642180822d0,
1298
# 1.311566449d0, 1.458302972d0, 1.610089467d0, 1.771001228d0,
1299
# 1.958842368d0, 1.980856139d0, 2.003769992d0, 2.027683661d0,
1300
# 2.052701483d0, 2.078930176d0, 2.106476183d0, 2.135441928d0,
1301
# 2.165921345d0, 2.197994437d0, 2.231721476d0, 2.323417121d0,
1302
# 2.425109533d0, 2.534671123d0, 2.648366651d0, 2.869451789d0,
1303
# 3.056555683d0, 3.197055532d0, 3.294656543d0, 3.399442397d0,
1304
# 3.438265245d0, 3.450036365d0, 3.449871888d0, 3.442815310d0,
1305
# 3.428262691d0, 3.411797242d0, 3.363961271d0, 3.295179282d0,
1306
# 3.287861332d0, 3.325219720d0, 3.470009593d0, 3.737173473d0,
1307
# 4.021846379d0, 4.241836705d0, 4.394931599d0, 4.501421364d0,
1308
# 4.577551451d0, 4.633363121d0, 4.675411596d0, 4.707771996d0,
1309
# 4.732593619d0, 4.751501682d0, 4.766245197d0, 4.778023555d0,
1310
# 4.787199047d0, 4.793849566d0, 4.798265576d0, 4.801110206d0,
1311
# 4.802711433d0, 4.767182955d0, 4.677220707d0, 4.558671105d0,
1312
# 4.423406368d0, 4.276574444d0, 4.119375511d0, 3.974083777d0,
1313
# 3.819733346d0, 3.671217773d0, 3.517370222d0, 3.344189966d0,
1314
# 3.116802834d0, 2.754066836d0, 1.809839704d0, 1.691727991d0,
1315
# 1.607871864d0, 1.544493194d0, 1.473234377d0, 1.390992309d0,
1316
# 1.292466985d0, 1.234537350d0, 1.167169430d0, 1.128863846d0,
1317
# 1.085963030d0, 1.037416400d0, 0.979461824d0, 0.952439442d0,
1318
# 0.922628624d0, 0.888386832d0, 0.846809865d0, 0.791742524d0,
1319
# 0.752347274d0, 0.742175952d0, 0.730422928d0, 0.716838887d0,
1320
# 0.697956938d0, 0.684261787d0, 0.672163545d0, 0.665924128d0,
1321
# 0.661453423d0, 0.652439216d0, 0.656998002d0, 0.677275168d0,
1322
# 0.699459467d0, 0.721833458d0, 0.738243489d0, 0.751895487d0,
1323
# 0.769698326d0, 0.779779847d0, 0.799512664d0, 0.816893187d0,
1324
# 0.842794338d0, 0.874519045d0, 0.916496501d0, 0.950264630d0,
1325
# 0.968728477d0, 1.105576625d0, 1.278360270d0, 1.391736231d0,
1326
# 1.542276157d0, 1.675869257d0, 1.795595314d0, 1.923303286d0,
1327
# 2.025352084d0, 2.157383377d0, 2.260417650d0, 2.361640519d0,
1328
# 2.479888170d0, 2.572642417d0, 2.688809606d0, 2.796999182d0,
1329
# 2.879392617d0, 3.848596335d0, 4.732837134d0, 5.525766484d0,
1330
# 6.243989288d0, 6.975830059d0, 7.607337162d0, 8.167782071d0,
1331
# 8.790374186d0, 9.422028844d0, 9.922844338d0, 10.426856521d0,
1332
# 10.933372206d0, 11.329434560d0, 11.735173309d0, 11.795158880d0,
1333
# 11.720854043d0, 11.455441492d0, 10.516061928d0, 9.017310079d0,
1334
# 7.000279246d0, 5.560369962d0, 3.715234200d0, 1.373399680d0,
1335
# -1.698051464d0, -5.394210410d0, -9.477907268d0,-13.677692505d0,
1336
# -18.256012300d0,-24.029072069d0,-29.992304980d0,-35.639341296d0,
1337
# -41.517060283d0,-46.430415213d0,-50.471411203d0,-53.247623713d0,
1338
# -55.184140125d0,-55.523807703d0,-55.749322940d0,-55.785869793d0,
1339
# -56.086101947d0,-56.289055861d0,-56.369321707d0,-56.555818993d0,
1340
# -56.632365789d0,-56.616832544d0,-56.709629986d0,-56.539655434d0,
1341
# -56.366567734d0,-56.272247546d0,-56.170833793d0,-55.984866548d0,
1342
# -55.776968735d0,-55.354292076d0,-55.136594171d0,-54.856130174d0,
1343
# -54.390381982d0,-50.585126291d0,-46.161103219d0,-41.187366158d0,
1344
# -36.554683053d0,-32.185245783d0,-28.097939638d0,-24.608441494d0,
1345
# -21.643195994d0,-18.984106588d0,-16.998672281d0,-15.273879570d0,
1346
# -13.373334219d0,-11.294832079d0, -9.584938331d0, -8.286389135d0,
1347
# -6.838755740d0, -5.409090746d0, -4.588339018d0, -3.705250749d0,
1348
# -2.476432078d0, -1.367536281d0, -0.414988031d0, 0.510347996d0,
1349
# 1.376925204d0, 2.276754999d0, 3.203065379d0, 3.629617497d0,
1350
# 3.821185286d0, 4.122178708d0, 4.804270233d0, 5.929961061d0,
1351
# 7.072311664d0, 7.458878393d0, 7.823110180d0/
1352
*
1353
DATA (xc(i),i=1,223)/100.0000d0,105.0000d0,110.0000d0,
1354
# 115.0000d0,120.0000d0,125.0000d0,130.0000d0,135.0000d0,
1355
# 140.0000d0,145.0000d0,150.0000d0,155.0000d0,156.0000d0,
1356
# 157.0000d0,158.0000d0,159.0000d0,159.1000d0,159.2000d0,
1357
# 159.3000d0,159.4000d0,159.5000d0,159.6000d0,159.7000d0,
1358
# 159.8000d0,159.9000d0,160.0000d0,160.2500d0,160.5000d0,
1359
# 160.7500d0,161.0000d0,161.5000d0,162.0000d0,162.5000d0,
1360
# 163.0000d0,164.0000d0,165.0000d0,166.0000d0,167.0000d0,
1361
# 168.0000d0,169.0000d0,170.0000d0,175.0000d0,178.5000d0,
1362
# 179.0000d0,179.5000d0,180.0000d0,180.5000d0,181.0000d0,
1363
# 181.5000d0,182.0000d0,182.5000d0,183.0000d0,184.0000d0,
1364
# 185.0000d0,186.0000d0,187.0000d0,188.0000d0,189.0000d0,
1365
# 190.0000d0,191.0000d0,192.0000d0,193.0000d0,194.0000d0,
1366
# 195.0000d0,196.0000d0,197.0000d0,198.0000d0,199.0000d0,
1367
# 200.0000d0,210.0000d0,220.0000d0,230.0000d0,240.0000d0,
1368
# 250.0000d0,260.0000d0,270.0000d0,280.0000d0,290.0000d0,
1369
# 300.0000d0,310.0000d0,320.0000d0,330.0000d0,340.0000d0,
1370
# 341.0000d0,342.0000d0,343.0000d0,344.0000d0,344.2000d0,
1371
# 344.4000d0,344.6000d0,344.8000d0,344.9000d0,345.0000d0,
1372
# 345.0500d0,345.1000d0,345.1500d0,345.2000d0,345.2200d0,
1373
# 345.2400d0,345.2600d0,345.2700d0,345.2720d0,345.2740d0,
1374
# 345.2760d0,345.2780d0,345.2790d0,345.2796d0,345.2798d0,
1375
# 345.2799d0,345.2804d0,345.2820d0,345.2900d0,345.3500d0,
1376
# 345.4500d0,345.5200d0,345.6000d0,345.7000d0,345.8000d0,
1377
# 346.0000d0,347.0000d0,348.0000d0,349.0000d0,350.0000d0,
1378
# 351.0000d0,352.0000d0,353.0000d0,354.0000d0,355.0000d0,
1379
# 356.0000d0,357.0000d0,358.0000d0,359.0000d0,360.0000d0,
1380
# 370.0000d0,380.0000d0,390.0000d0,400.0000d0,410.0000d0,
1381
# 420.0000d0,430.0000d0,440.0000d0,450.0000d0,460.0000d0,
1382
# 470.0000d0,480.0000d0,490.0000d0,500.0000d0,510.0000d0,
1383
# 520.0000d0,530.0000d0,540.0000d0,550.0000d0,560.0000d0,
1384
# 565.0000d0,570.0000d0,575.0000d0,580.0000d0,585.0000d0,
1385
# 590.0000d0,595.0000d0,600.0000d0,605.0000d0,610.0000d0,
1386
# 615.0000d0,620.0000d0,625.0000d0,630.0000d0,635.0000d0,
1387
# 640.0000d0,641.0000d0,642.0000d0,643.0000d0,644.0000d0,
1388
# 645.0000d0,646.0000d0,647.0000d0,648.0000d0,649.0000d0,
1389
# 650.0000d0,651.0000d0,652.0000d0,653.0000d0,654.0000d0,
1390
# 655.0000d0,656.0000d0,657.0000d0,658.0000d0,659.0000d0,
1391
# 660.0000d0,670.0000d0,680.0000d0,690.0000d0,700.0000d0,
1392
# 710.0000d0,720.0000d0,730.0000d0,740.0000d0,750.0000d0,
1393
# 760.0000d0,770.0000d0,780.0000d0,790.0000d0,800.0000d0,
1394
# 810.0000d0,820.0000d0,830.0000d0,840.0000d0,850.0000d0,
1395
# 860.0000d0,870.0000d0,880.0000d0,890.0000d0,900.0000d0,
1396
# 910.0000d0,920.0000d0,930.0000d0,940.0000d0,950.0000d0,
1397
# 960.0000d0,970.0000d0,980.0000d0,990.0000d0,1000.0000d0/
1398
*
1399
DATA (yc(i),i=1,223)/-2.969458664d0,-2.755568498d0,-2.518689112d0,
1400
# -2.256840599d0, -1.967703151d0, -1.648299088d0, -1.294794285d0,
1401
# -0.902525232d0, -0.463541061d0, 0.032517848d0, 0.603841958d0,
1402
# 1.273887882d0, 1.420789352d0, 1.572747980d0, 1.733843353d0,
1403
# 1.921879203d0, 1.943912705d0, 1.966846336d0, 1.990779880d0,
1404
# 2.015817728d0, 2.042066621d0, 2.069633010d0, 2.098619322d0,
1405
# 2.129119489d0, 2.161213532d0, 2.194961744d0, 2.286711412d0,
1406
# 2.388459445d0, 2.498078035d0, 2.611830784d0, 2.833024715d0,
1407
# 3.020227220d0, 3.160820150d0, 3.258502948d0, 3.363381747d0,
1408
# 3.402305539d0, 3.414208425d0, 3.414049050d0, 3.407053730d0,
1409
# 3.392629803d0, 3.376195597d0, 3.328499111d0, 3.268651902d0,
1410
# 3.259860594d0, 3.253486002d0, 3.252546948d0, 3.262154622d0,
1411
# 3.289924633d0, 3.345085115d0, 3.434732330d0, 3.557729513d0,
1412
# 3.701897111d0, 3.986569984d0, 4.206557658d0, 4.359682402d0,
1413
# 4.466176116d0, 4.542300703d0, 4.598184529d0, 4.640273713d0,
1414
# 4.672604523d0, 4.697438255d0, 4.716361927d0, 4.731103154d0,
1415
# 4.742925539d0, 4.752092918d0, 4.758670785d0, 4.763130070d0,
1416
# 4.766074843d0, 4.767748196d0, 4.732269356d0, 4.642413616d0,
1417
# 4.524181528d0, 4.389510854d0, 4.242822112d0, 4.086217524d0,
1418
# 3.941680893d0, 3.789691168d0, 3.644241595d0, 3.494351920d0,
1419
# 3.331867639d0, 3.122929519d0, 2.816374576d0, 2.186038507d0,
1420
# 2.068538445d0, 1.926511260d0, 1.746359564d0, 1.494945494d0,
1421
# 1.428923128d0, 1.354291504d0, 1.267892764d0, 1.163289817d0,
1422
# 1.100587592d0, 1.027248191d0, 0.984553881d0, 0.936444199d0,
1423
# 0.879742362d0, 0.808975906d0, 0.774145495d0, 0.731724792d0,
1424
# 0.675557165d0, 0.635406715d0, 0.625004410d0, 0.613030568d0,
1425
# 0.599164524d0, 0.579800865d0, 0.565857087d0, 0.553503488d0,
1426
# 0.547057831d0, 0.542411799d0, 0.533743538d0, 0.538671942d0,
1427
# 0.542814415d0, 0.562187598d0, 0.575771207d0, 0.586370628d0,
1428
# 0.605546923d0, 0.622661850d0, 0.639021282d0, 0.683852301d0,
1429
# 0.846885709d0, 0.995547199d0, 1.139428677d0, 1.270752033d0,
1430
# 1.403862461d0, 1.540060220d0, 1.676615197d0, 1.791154110d0,
1431
# 1.909977528d0, 2.023870438d0, 2.138393119d0, 2.233199625d0,
1432
# 2.358706648d0, 2.446271802d0, 3.459751114d0, 4.411864477d0,
1433
# 5.220622691d0, 5.988967245d0, 6.706065217d0, 7.359912016d0,
1434
# 7.993290467d0, 8.675262256d0, 9.327550412d0, 9.917129519d0,
1435
# 10.373037136d0, 11.009210530d0, 11.516244354d0, 11.985967329d0,
1436
# 12.198319618d0, 12.357011718d0, 12.311558072d0, 11.754465713d0,
1437
# 10.617092695d0, 9.043395280d0, 8.002055503d0, 6.601203443d0,
1438
# 4.587232698d0, 1.906399387d0, -1.281177539d0, -5.052281649d0,
1439
# -9.213421228d0,-13.868726558d0,-19.489194521d0,-25.988103189d0,
1440
# -32.734525726d0,-39.402198052d0,-45.903352541d0,-51.189352305d0,
1441
# -55.322071720d0,-58.186905136d0,-58.653812374d0,-59.048160142d0,
1442
# -59.279596622d0,-59.992120796d0,-60.271154555d0,-60.549383981d0,
1443
# -60.729563975d0,-60.954540486d0,-61.084961955d0,-61.164811836d0,
1444
# -61.186911575d0,-61.109847454d0,-61.108201824d0,-60.992397754d0,
1445
# -60.768938455d0,-60.645603901d0,-60.224362756d0,-59.860911894d0,
1446
# -59.548720101d0,-59.301307840d0,-55.034929471d0,-49.745997496d0,
1447
# -44.181440809d0,-39.061582838d0,-34.158470096d0,-29.603188102d0,
1448
# -25.799341483d0,-22.634533844d0,-19.804942358d0,-17.654326276d0,
1449
# -15.913630982d0,-13.830699918d0,-11.633319645d0,-10.006868842d0,
1450
# -8.532937676d0, -7.132846889d0, -5.627566626d0, -4.758444504d0,
1451
# -3.887401469d0, -2.729580751d0, -1.512994295d0, -0.550651909d0,
1452
# 0.355022046d0, 1.248689020d0, 2.123963438d0, 3.014801413d0,
1453
# 3.484406765d0, 3.696913629d0, 4.008426964d0, 4.651255618d0,
1454
# 5.840647264d0, 6.935940541d0, 7.364646872d0, 7.670905366d0/
1455
*
1456
DATA (xu(i),i=1,225)/100.0000d0,105.0000d0,110.0000d0,
1457
# 115.0000d0,120.0000d0,125.0000d0,130.0000d0,135.0000d0,
1458
# 140.0000d0,145.0000d0,150.0000d0,155.0000d0,156.0000d0,
1459
# 157.0000d0,158.0000d0,159.0000d0,159.1000d0,159.2000d0,
1460
# 159.3000d0,159.4000d0,159.5000d0,159.6000d0,159.7000d0,
1461
# 159.8000d0,159.9000d0,160.0000d0,160.2500d0,160.5000d0,
1462
# 160.7500d0,161.0000d0,161.5000d0,162.0000d0,162.5000d0,
1463
# 163.0000d0,164.0000d0,165.0000d0,166.0000d0,167.0000d0,
1464
# 168.0000d0,169.0000d0,170.0000d0,175.0000d0,178.5000d0,
1465
# 179.0000d0,179.5000d0,180.0000d0,180.5000d0,181.0000d0,
1466
# 181.5000d0,182.0000d0,182.5000d0,183.0000d0,184.0000d0,
1467
# 185.0000d0,186.0000d0,187.0000d0,188.0000d0,189.0000d0,
1468
# 190.0000d0,191.0000d0,192.0000d0,193.0000d0,194.0000d0,
1469
# 195.0000d0,196.0000d0,197.0000d0,198.0000d0,199.0000d0,
1470
# 200.0000d0,210.0000d0,220.0000d0,230.0000d0,240.0000d0,
1471
# 250.0000d0,260.0000d0,270.0000d0,280.0000d0,290.0000d0,
1472
# 300.0000d0,310.0000d0,320.0000d0,330.0000d0,340.0000d0,
1473
# 341.0000d0,342.0000d0,343.0000d0,344.0000d0,345.0000d0,
1474
# 345.5000d0,345.8000d0,346.0000d0,347.0000d0,347.2000d0,
1475
# 347.4000d0,347.6000d0,347.8000d0,348.0000d0,348.1000d0,
1476
# 348.2000d0,348.2500d0,348.3000d0,348.3200d0,348.3400d0,
1477
# 348.3600d0,348.3800d0,348.4000d0,348.4100d0,348.4200d0,
1478
# 348.4300d0,348.4320d0,348.4340d0,348.4360d0,348.4380d0,
1479
# 348.4390d0,348.4396d0,348.4398d0,348.4399d0,348.4404d0,
1480
# 348.4420d0,348.4500d0,348.5000d0,348.7000d0,349.0000d0,
1481
# 349.4000d0,350.0000d0,351.0000d0,352.0000d0,353.0000d0,
1482
# 354.0000d0,355.0000d0,356.0000d0,357.0000d0,358.0000d0,
1483
# 359.0000d0,360.0000d0,370.0000d0,380.0000d0,390.0000d0,
1484
# 400.0000d0,410.0000d0,420.0000d0,430.0000d0,440.0000d0,
1485
# 450.0000d0,460.0000d0,470.0000d0,480.0000d0,490.0000d0,
1486
# 500.0000d0,510.0000d0,520.0000d0,530.0000d0,540.0000d0,
1487
# 550.0000d0,560.0000d0,565.0000d0,570.0000d0,575.0000d0,
1488
# 580.0000d0,585.0000d0,590.0000d0,595.0000d0,600.0000d0,
1489
# 605.0000d0,610.0000d0,615.0000d0,620.0000d0,625.0000d0,
1490
# 630.0000d0,635.0000d0,640.0000d0,641.0000d0,642.0000d0,
1491
# 643.0000d0,644.0000d0,645.0000d0,646.0000d0,647.0000d0,
1492
# 648.0000d0,649.0000d0,650.0000d0,651.0000d0,652.0000d0,
1493
# 653.0000d0,654.0000d0,655.0000d0,656.0000d0,657.0000d0,
1494
# 658.0000d0,659.0000d0,660.0000d0,670.0000d0,680.0000d0,
1495
# 690.0000d0,700.0000d0,710.0000d0,720.0000d0,730.0000d0,
1496
# 740.0000d0,750.0000d0,760.0000d0,770.0000d0,780.0000d0,
1497
# 790.0000d0,800.0000d0,810.0000d0,820.0000d0,830.0000d0,
1498
# 840.0000d0,850.0000d0,860.0000d0,870.0000d0,880.0000d0,
1499
# 890.0000d0,900.0000d0,910.0000d0,920.0000d0,930.0000d0,
1500
# 940.0000d0,950.0000d0,960.0000d0,970.0000d0,980.0000d0,
1501
# 990.0000d0,1000.0000d0/
1502
*
1503
DATA (yu(i),i=1,225)/-3.011799533d0,-2.797658814d0,-2.560511342d0,
1504
# -2.298374579d0, -2.008927458d0, -1.689194746d0, -1.335334782d0,
1505
# -0.942678003d0, -0.503267763d0, -0.006741983d0, 0.565124142d0,
1506
# 1.235840116d0, 1.382905273d0, 1.535044791d0, 1.696328604d0,
1507
# 1.884567714d0, 1.906622347d0, 1.929577239d0, 1.953532153d0,
1508
# 1.978591436d0, 2.004861792d0, 2.032449634d0, 2.061457355d0,
1509
# 2.091978868d0, 2.124094156d0, 2.157863499d0, 2.249665540d0,
1510
# 2.351465763d0, 2.461136573d0, 2.574941494d0, 2.796240093d0,
1511
# 2.983545301d0, 3.124228166d0, 3.221982342d0, 3.326972382d0,
1512
# 3.365993805d0, 3.378023193d0, 3.377887062d0, 3.370931083d0,
1513
# 3.356634680d0, 3.340265619d0, 3.292682936d0, 3.232961811d0,
1514
# 3.224194020d0, 3.217828260d0, 3.216889093d0, 3.226496714d0,
1515
# 3.254271720d0, 3.309440914d0, 3.399097855d0, 3.522103315d0,
1516
# 3.666276367d0, 3.950949463d0, 4.170926862d0, 4.324070476d0,
1517
# 4.430588837d0, 4.506706461d0, 4.562622634d0, 4.604788444d0,
1518
# 4.637112788d0, 4.661928180d0, 4.680880990d0, 4.695625131d0,
1519
# 4.707459275d0, 4.716670603d0, 4.723178694d0, 4.727608289d0,
1520
# 4.730650318d0, 4.732442748d0, 4.697030978d0, 4.607168713d0,
1521
# 4.489397175d0, 4.355888392d0, 4.208987612d0, 4.052707904d0,
1522
# 3.908019896d0, 3.759968739d0, 3.615348691d0, 3.469932507d0,
1523
# 3.315332978d0, 3.122286809d0, 2.855004727d0, 2.377722123d0,
1524
# 2.300655190d0, 2.213908305d0, 2.115631830d0, 2.000162289d0,
1525
# 1.860840425d0, 1.779119559d0, 1.724912832d0, 1.685992661d0,
1526
# 1.446048759d0, 1.384730617d0, 1.316264190d0, 1.238304389d0,
1527
# 1.147571101d0, 1.036738128d0, 0.969143795d0, 0.888640045d0,
1528
# 0.841031364d0, 0.785526212d0, 0.759880178d0, 0.731970096d0,
1529
# 0.701136254d0, 0.665654946d0, 0.622492994d0, 0.596679361d0,
1530
# 0.565315972d0, 0.524528227d0, 0.513928834d0, 0.501826683d0,
1531
# 0.487735928d0, 0.468019598d0, 0.454013349d0, 0.441470198d0,
1532
# 0.434889761d0, 0.430148561d0, 0.421464835d0, 0.429200773d0,
1533
# 0.429255984d0, 0.437861500d0, 0.480919141d0, 0.544499260d0,
1534
# 0.606035270d0, 0.710025546d0, 0.864539937d0, 1.009890014d0,
1535
# 1.156589498d0, 1.268264088d0, 1.408793136d0, 1.528293889d0,
1536
# 1.642389803d0, 1.771929787d0, 1.875426734d0, 1.979445634d0,
1537
# 3.092059170d0, 4.047318673d0, 4.937427406d0, 5.684883866d0,
1538
# 6.473364254d0, 7.141842873d0, 7.823212681d0, 8.484967861d0,
1539
# 9.214907909d0, 9.784981877d0, 10.383431544d0, 11.054527450d0,
1540
# 11.588357911d0, 12.165403348d0, 12.507103232d0, 12.848157153d0,
1541
# 12.975158246d0, 12.770167330d0, 11.991429721d0, 11.123821384d0,
1542
# 10.241165300d0, 9.157568585d0, 7.611314332d0, 5.499567320d0,
1543
# 2.603705830d0, -0.579201122d0, -4.348429690d0, -8.823492981d0,
1544
# -14.416229606d0,-20.848895478d0,-28.011934618d0,-36.077657061d0,
1545
# -43.726529167d0,-50.682419374d0,-56.593917330d0,-60.848713817d0,
1546
# -61.464423821d0,-62.283661600d0,-62.930084228d0,-63.549395305d0,
1547
# -64.133867221d0,-64.465495513d0,-64.909336911d0,-65.479678788d0,
1548
# -65.703953492d0,-66.113782243d0,-66.222507344d0,-66.278786408d0,
1549
# -66.393237110d0,-66.366867383d0,-66.202199439d0,-66.029041364d0,
1550
# -65.774112798d0,-65.485847916d0,-65.247943137d0,-64.801055924d0,
1551
# -60.101647018d0,-54.084245796d0,-47.571709079d0,-41.661629810d0,
1552
# -36.168693788d0,-31.212363345d0,-27.000704767d0,-23.605327875d0,
1553
# -20.580041179d0,-18.274370612d0,-16.434309906d0,-14.327218251d0,
1554
# -12.021979554d0,-10.331796926d0, -8.883031936d0, -7.348435651d0,
1555
# -5.780362152d0, -4.953756827d0, -4.070387325d0, -2.840184650d0,
1556
# -1.678475825d0, -0.679391004d0, 0.223829200d0, 1.086204496d0,
1557
# 2.026633865d0, 2.923942775d0, 3.306978920d0, 3.579966278d0,
1558
# 3.863101723d0, 4.490999475d0, 5.727923150d0, 6.829991246d0,
1559
# 7.212458092d0, 7.551889308d0/
1560
*
1561
IF(top.eq.-1) THEN
1562
n= 226
1563
DO l=1,226
1564
x(l)= xl(l)
1565
y(l)= yl(l)
1566
c write(35,*)l,x(l),y(l)
1567
ENDDO
1568
ELSEIF(top.eq.0) THEN
1569
n= 223
1570
DO l=1,223
1571
x(l)= xc(l)
1572
y(l)= yc(l)
1573
c write(36,*)l,x(l),y(l)
1574
ENDDO
1575
ELSEIF(top.eq.1) THEN
1576
n= 225
1577
DO l=1,225
1578
x(l)= xu(l)
1579
y(l)= yu(l)
1580
c write(37,*)l,x(l),y(l)
1581
ENDDO
1582
ENDIF
1583
*
1584
*.....First check if u is in the same interval found on the
1585
* last call to Seval.............................................
1586
*
1587
IF ( (i.lt.1) .OR. (i.ge.n) ) i=1
1588
IF ( (u.lt.x(i)) .OR. (u.ge.x(i+1)) ) THEN
1589
i=1
1590
*
1591
* binary search
1592
*
1593
j= n+1
1594
2468 k= (i+j)/2
1595
IF (u.lt.x(k)) THEN
1596
j= k
1597
ELSE
1598
i= k
1599
ENDIF
1600
IF (j.le.i+1) GOTO 2469
1601
GOTO 2468
1602
2469 CONTINUE
1603
ENDIF
1604
*
1605
dx= u-x(i)
1606
*
1607
* evaluate the spline
1608
*
1609
f= y(i)+dx*(b(i)+dx*(c(i)+dx*d(i)))
1610
fp= b(i)+dx*(TWO*c(i) + dx*THREE*d(i))
1611
fpp= TWO*c(i) + dx*SIX*d(i)
1612
fppp= SIX*d(i)
1613
*
1614
RETURN
1615
*
1616
END
1617
1618
*****************************************************************
1619
C---electroweak corrections to H --> gg in SM4
1620
C m_t' = m_b' + 50*(1+1/5*log(M_H/(115 GeV)) GeV
1621
*****************************************************************
1622
1623
double precision function glgl_elw4(im,mh)
1624
*
1625
IMPLICIT NONE
1626
*
1627
INTEGER i,top,gdim,im
1628
REAL*8 u,ui,value
1629
REAL*8 vp,vpp,vppp
1630
REAL*8 mtop,mh
1631
REAL*8 x,x0,x1,x2,y0,y1,y2,a0,a1,a2
1632
REAL*8 value0,value1,value2
1633
REAL*8 bl(94),cl(94),dl(94)
1634
REAL*8 bc(120),cc(120),dc(120)
1635
*
1636
* u value of M_H at which the spline is to be evaluated
1637
* top= -1,0,1 lower, central, upper value for m_top
1638
*
1639
u = mh
1640
if(im.eq.1)then
1641
gdim= 94
1642
CALL F4MMsplineSingle1(bl,cl,dl,gdim)
1643
CALL S4eval3Single1(u,bl,cl,dl,gdim,value,vp,vpp,vppp)
1644
else
1645
gdim= 120
1646
CALL F4MMsplineSingle2(bc,cc,dc,gdim)
1647
CALL S4eval3Single2(u,bc,cc,dc,gdim,value,vp,vpp,vppp)
1648
endif
1649
glgl_elw4=value/100.d0
1650
RETURN
1651
END
1652
*
1653
*-----------------------------------------------------------------------
1654
*
1655
c CONTAINS
1656
*
1657
SUBROUTINE F4MMsplineSingle1(b,c,d,gdim)
1658
*
1659
* ---------------------------------------------------------------------------
1660
*
1661
* PURPOSE - Compute the coefficients b,c,d for a cubic interpolating spline
1662
* so that the interpolated value is given by
1663
* s(x) = y(k) + b(k)*(x-x(k)) + c(k)*(x-x(k))**2 + d(k)*(x-x(k))**3
1664
* when x(k) <= x <= x(k+1)
1665
* The end conditions match the third derivatives of the interpolated curve to
1666
* the third derivatives of the unique polynomials thru the first four and
1667
* last four points.
1668
* Use Seval or Seval3 to evaluate the spline.
1669
*
1670
INTEGER k,n,i,top,gdim,l
1671
*
1672
REAL*8 xc(94),yc(94)
1673
REAL*8 x(94),y(94)
1674
*
1675
REAL*8 b(gdim)
1676
* linear coeff
1677
*
1678
REAL*8 c(gdim)
1679
* quadratic coeff.
1680
*
1681
REAL*8 d(gdim)
1682
* cubic coeff.
1683
*
1684
REAL*8 t
1685
REAL*8 ZERO, TWO, THREE
1686
PARAMETER(ZERO=0.0, TWO=2.0, THREE=3.0)
1687
*
1688
* The grid
1689
*
1690
DATA (xc(i),i=1,94)/
1691
.100.000000000, 110.000000000, 120.000000000, 130.000000000,
1692
.140.000000000, 145.000000000, 150.000000000, 151.000000000,
1693
.152.000000000, 153.000000000, 154.000000000, 155.000000000,
1694
.156.000000000, 157.000000000, 158.000000000, 159.000000000,
1695
.160.000000000, 161.000000000, 162.000000000, 163.000000000,
1696
.164.000000000, 165.000000000, 166.000000000, 167.000000000,
1697
.168.000000000, 169.000000000, 170.000000000, 171.000000000,
1698
.172.000000000, 173.000000000, 174.000000000, 175.000000000,
1699
.176.000000000, 177.000000000, 178.000000000, 179.000000000,
1700
.180.000000000, 181.000000000, 182.000000000, 183.000000000,
1701
.184.000000000, 185.000000000, 186.000000000, 187.000000000,
1702
.188.000000000, 189.000000000, 190.000000000, 195.000000000,
1703
.200.000000000, 210.000000000, 220.000000000, 230.000000000,
1704
.240.000000000, 250.000000000, 260.000000000, 270.000000000,
1705
.280.000000000, 290.000000000, 300.000000000, 310.000000000,
1706
.325.000000000, 330.000000000, 335.000000000, 337.000000000,
1707
.339.000000000, 341.000000000, 342.000000000, 343.000000000,
1708
.344.000000000, 345.000000000, 346.000000000, 347.000000000,
1709
.348.000000000, 349.000000000, 350.000000000, 355.000000000,
1710
.360.000000000, 365.000000000, 370.000000000, 380.000000000,
1711
.390.000000000, 395.000000000, 400.000000000, 420.000000000,
1712
.430.000000000, 440.000000000, 450.000000000, 460.000000000,
1713
.470.000000000, 480.000000000, 490.000000000, 500.000000000,
1714
.550.000000000, 600.000000000/
1715
*
1716
DATA (yc(i),i=1,94)/
1717
. 7.081784858, 7.005040906, 6.909638450, 6.767663371,
1718
. 6.553481059, 6.394170663, 6.161453100, 6.100169493,
1719
. 6.031159987, 5.952322286, 5.860793898, 5.752611211,
1720
. 5.622148850, 5.462299560, 5.267015468, 5.045635485,
1721
. 4.867448493, 4.858996855, 4.957159679, 4.991788742,
1722
. 4.949709854, 4.869086008, 4.773624254, 4.672977010,
1723
. 4.574106306, 4.473957550, 4.375216682, 4.277835718,
1724
. 4.180106440, 4.082355068, 3.981303187, 3.876580238,
1725
. 3.765612293, 3.650290738, 3.517994735, 3.374532294,
1726
. 3.224190879, 3.081402574, 3.003176187, 3.023510652,
1727
. 3.058333994, 3.061527486, 3.032524574, 2.980103496,
1728
. 2.916221523, 2.850606479, 2.785863031, 2.465569133,
1729
. 2.204423840, 1.758167155, 1.356606770, 1.009560117,
1730
. 0.701093096, 0.393933499, 0.079112288, -0.206013527,
1731
.-0.518586518, -0.811157249, -1.129144131, -1.439334674,
1732
.-1.957798896, -2.212619484, -2.441510248, -2.573014363,
1733
.-2.728219541, -2.899530309, -2.993431692, -3.158243996,
1734
.-3.351666868, -3.920928883, -4.098697186, -4.132701122,
1735
.-4.144213562, -4.088523452, -4.049143398, -3.913706462,
1736
.-3.749072205, -3.669404489, -3.618383929, -3.575928277,
1737
.-3.688589851, -3.744277795, -3.841702830, -4.446704386,
1738
.-4.761423299, -5.287974189, -5.744301071, -6.211157195,
1739
.-6.778672157, -7.445798763, -8.061705151, -8.714706599,
1740
.-12.512297879, -16.999567032/
1741
*
1742
n= 94
1743
DO l=1,n
1744
x(l)= xc(l)
1745
y(l)= yc(l)
1746
ENDDO
1747
1748
*.....Set up tridiagonal system.........................................
1749
* b=diagonal, d=offdiagonal, c=right-hand side
1750
*
1751
d(1)= x(2)-x(1)
1752
c(2)= (y(2)-y(1))/d(1)
1753
DO k= 2,n-1
1754
d(k)= x(k+1)-x(k)
1755
b(k)= TWO*(d(k-1)+d(k))
1756
c(k+1)= (y(k+1)-y(k))/d(k)
1757
c(k)= c(k+1)-c(k)
1758
END DO
1759
*
1760
*.....End conditions. third derivatives at x(1) and x(n) obtained
1761
* from divided differences.......................................
1762
*
1763
b(1)= -d(1)
1764
b(n)= -d(n-1)
1765
c(1)= ZERO
1766
c(n)= ZERO
1767
IF (n.gt.3) THEN
1768
c(1)= c(3)/(x(4)-x(2))-c(2)/(x(3)-x(1))
1769
c(n)= c(n-1)/(x(n)-x(n-2))-c(n-2)/(x(n-1)-x(n-3))
1770
c(1)= c(1)*d(1)*d(1)/(x(4)-x(1))
1771
c(n)= -c(n)*d(n-1)*d(n-1)/(x(n)-x(n-3))
1772
END IF
1773
*
1774
DO k=2,n ! forward elimination
1775
t= d(k-1)/b(k-1)
1776
b(k)= b(k)-t*d(k-1)
1777
c(k)= c(k)-t*c(k-1)
1778
END DO
1779
*
1780
c(n)= c(n)/b(n)
1781
*
1782
* back substitution ( makes c the sigma of text)
1783
*
1784
DO k=n-1,1,-1
1785
c(k)= (c(k)-d(k)*c(k+1))/b(k)
1786
END DO
1787
*
1788
*.....Compute polynomial coefficients...................................
1789
*
1790
b(n)= (y(n)-y(n-1))/d(n-1)+d(n-1)*(c(n-1)+c(n)+c(n))
1791
DO k=1,n-1
1792
b(k)= (y(k+1)-y(k))/d(k)-d(k)*(c(k+1)+c(k)+c(k))
1793
d(k)= (c(k+1)-c(k))/d(k)
1794
c(k)= THREE*c(k)
1795
END DO
1796
c(n)= THREE*c(n)
1797
d(n)= d(n-1)
1798
*
1799
RETURN
1800
*
1801
END
1802
*
1803
*------------------------------------------------------------------------
1804
*
1805
SUBROUTINE S4eval3Single1(u,b,c,d,gdim,f,fp,fpp,fppp)
1806
*
1807
* ---------------------------------------------------------------------------
1808
*
1809
* PURPOSE - Evaluate the cubic spline function
1810
* Seval=y(i)+b(i)*(u-x(i))+c(i)*(u-x(i))**2+d(i)*(u-x(i))**3
1811
* where x(i) <= u < x(i+1)
1812
*
1813
* NOTES- if u<x(1), i=1 is used;if u>x(n), i=n is used
1814
*
1815
REAL*8 u
1816
* abscissa at which the spline is to be evaluated
1817
*
1818
INTEGER j,k,n,l,top,gdim
1819
*
1820
REAL*8 xc(94),yc(94)
1821
REAL*8 x(94),y(94)
1822
REAL*8 b(gdim),c(gdim),d(gdim)
1823
* linear,quadratic,cubic coeff
1824
*
1825
REAL*8 f,fp,fpp,fppp
1826
* function, 1st,2nd,3rd deriv
1827
*
1828
INTEGER i
1829
REAL*8 dx
1830
REAL*8 ZERO, TWO, THREE
1831
PARAMETER(ZERO=0.0, TWO=2.0, THREE=3.0)
1832
SAVE i
1833
DATA i/1/
1834
*
1835
* The grid
1836
*
1837
DATA (xc(l),l=1,94)/
1838
.100.000000000, 110.000000000, 120.000000000, 130.000000000,
1839
.140.000000000, 145.000000000, 150.000000000, 151.000000000,
1840
.152.000000000, 153.000000000, 154.000000000, 155.000000000,
1841
.156.000000000, 157.000000000, 158.000000000, 159.000000000,
1842
.160.000000000, 161.000000000, 162.000000000, 163.000000000,
1843
.164.000000000, 165.000000000, 166.000000000, 167.000000000,
1844
.168.000000000, 169.000000000, 170.000000000, 171.000000000,
1845
.172.000000000, 173.000000000, 174.000000000, 175.000000000,
1846
.176.000000000, 177.000000000, 178.000000000, 179.000000000,
1847
.180.000000000, 181.000000000, 182.000000000, 183.000000000,
1848
.184.000000000, 185.000000000, 186.000000000, 187.000000000,
1849
.188.000000000, 189.000000000, 190.000000000, 195.000000000,
1850
.200.000000000, 210.000000000, 220.000000000, 230.000000000,
1851
.240.000000000, 250.000000000, 260.000000000, 270.000000000,
1852
.280.000000000, 290.000000000, 300.000000000, 310.000000000,
1853
.325.000000000, 330.000000000, 335.000000000, 337.000000000,
1854
.339.000000000, 341.000000000, 342.000000000, 343.000000000,
1855
.344.000000000, 345.000000000, 346.000000000, 347.000000000,
1856
.348.000000000, 349.000000000, 350.000000000, 355.000000000,
1857
.360.000000000, 365.000000000, 370.000000000, 380.000000000,
1858
.390.000000000, 395.000000000, 400.000000000, 420.000000000,
1859
.430.000000000, 440.000000000, 450.000000000, 460.000000000,
1860
.470.000000000, 480.000000000, 490.000000000, 500.000000000,
1861
.550.000000000, 600.000000000/
1862
*
1863
DATA (yc(l),l=1,94)/
1864
. 7.081784858, 7.005040906, 6.909638450, 6.767663371,
1865
. 6.553481059, 6.394170663, 6.161453100, 6.100169493,
1866
. 6.031159987, 5.952322286, 5.860793898, 5.752611211,
1867
. 5.622148850, 5.462299560, 5.267015468, 5.045635485,
1868
. 4.867448493, 4.858996855, 4.957159679, 4.991788742,
1869
. 4.949709854, 4.869086008, 4.773624254, 4.672977010,
1870
. 4.574106306, 4.473957550, 4.375216682, 4.277835718,
1871
. 4.180106440, 4.082355068, 3.981303187, 3.876580238,
1872
. 3.765612293, 3.650290738, 3.517994735, 3.374532294,
1873
. 3.224190879, 3.081402574, 3.003176187, 3.023510652,
1874
. 3.058333994, 3.061527486, 3.032524574, 2.980103496,
1875
. 2.916221523, 2.850606479, 2.785863031, 2.465569133,
1876
. 2.204423840, 1.758167155, 1.356606770, 1.009560117,
1877
. 0.701093096, 0.393933499, 0.079112288, -0.206013527,
1878
.-0.518586518, -0.811157249, -1.129144131, -1.439334674,
1879
.-1.957798896, -2.212619484, -2.441510248, -2.573014363,
1880
.-2.728219541, -2.899530309, -2.993431692, -3.158243996,
1881
.-3.351666868, -3.920928883, -4.098697186, -4.132701122,
1882
.-4.144213562, -4.088523452, -4.049143398, -3.913706462,
1883
.-3.749072205, -3.669404489, -3.618383929, -3.575928277,
1884
.-3.688589851, -3.744277795, -3.841702830, -4.446704386,
1885
.-4.761423299, -5.287974189, -5.744301071, -6.211157195,
1886
.-6.778672157, -7.445798763, -8.061705151, -8.714706599,
1887
.-12.512297879, -16.999567032/
1888
*
1889
n= 94
1890
DO l=1,n
1891
x(l)= xc(l)
1892
y(l)= yc(l)
1893
ENDDO
1894
*
1895
*.....First check if u is in the same interval found on the
1896
* last call to Seval.............................................
1897
*
1898
IF ( (i.lt.1) .OR. (i.ge.n) ) i=1
1899
IF ( (u.lt.x(i)) .OR. (u.ge.x(i+1)) ) THEN
1900
i=1
1901
*
1902
* binary search
1903
*
1904
j= n+1
1905
2468 k= (i+j)/2
1906
IF (u.lt.x(k)) THEN
1907
j= k
1908
ELSE
1909
i= k
1910
ENDIF
1911
IF (j.le.i+1) GOTO 2469
1912
GOTO 2468
1913
2469 CONTINUE
1914
ENDIF
1915
*
1916
dx= u-x(i)
1917
*
1918
* evaluate the spline
1919
*
1920
f= y(i)+dx*(b(i)+dx*(c(i)+dx*d(i)))
1921
fp= b(i)+dx*(TWO*c(i) + dx*THREE*d(i))
1922
fpp= TWO*c(i) + dx*SIX*d(i)
1923
fppp= SIX*d(i)
1924
*
1925
RETURN
1926
*
1927
END
1928
1929
SUBROUTINE F4MMsplineSingle2(b,c,d,gdim)
1930
*
1931
* ---------------------------------------------------------------------------
1932
*
1933
* PURPOSE - Compute the coefficients b,c,d for a cubic interpolating spline
1934
* so that the interpolated value is given by
1935
* s(x) = y(k) + b(k)*(x-x(k)) + c(k)*(x-x(k))**2 + d(k)*(x-x(k))**3
1936
* when x(k) <= x <= x(k+1)
1937
* The end conditions match the third derivatives of the interpolated curve to
1938
* the third derivatives of the unique polynomials thru the first four and
1939
* last four points.
1940
* Use Seval or Seval3 to evaluate the spline.
1941
*
1942
INTEGER k,n,i,top,gdim,l
1943
*
1944
REAL*8 xc(120),yc(120)
1945
REAL*8 x(120),y(120)
1946
*
1947
REAL*8 b(gdim)
1948
* linear coeff
1949
*
1950
REAL*8 c(gdim)
1951
* quadratic coeff.
1952
*
1953
REAL*8 d(gdim)
1954
* cubic coeff.
1955
*
1956
REAL*8 t
1957
REAL*8 ZERO, TWO, THREE
1958
PARAMETER(ZERO=0.0, TWO=2.0, THREE=3.0)
1959
*
1960
* The grid
1961
*
1962
DATA (xc(i),i=1,120)/
1963
.100.000000000, 110.000000000, 120.000000000, 130.000000000,
1964
.140.000000000, 145.000000000, 150.000000000, 151.000000000,
1965
.152.000000000, 153.000000000, 154.000000000, 155.000000000,
1966
.156.000000000, 157.000000000, 158.000000000, 159.000000000,
1967
.160.000000000, 161.000000000, 162.000000000, 163.000000000,
1968
.164.000000000, 165.000000000, 166.000000000, 167.000000000,
1969
.168.000000000, 169.000000000, 170.000000000, 171.000000000,
1970
.172.000000000, 173.000000000, 174.000000000, 175.000000000,
1971
.176.000000000, 177.000000000, 178.000000000, 179.000000000,
1972
.180.000000000, 181.000000000, 182.000000000, 183.000000000,
1973
.184.000000000, 185.000000000, 186.000000000, 187.000000000,
1974
.188.000000000, 189.000000000, 190.000000000, 195.000000000,
1975
.200.000000000, 210.000000000, 220.000000000, 230.000000000,
1976
.240.000000000, 250.000000000, 260.000000000, 270.000000000,
1977
.280.000000000, 290.000000000, 300.000000000, 310.000000000,
1978
.320.000000000, 325.000000000, 330.000000000, 335.000000000,
1979
.336.000000000, 337.000000000, 338.000000000, 339.000000000,
1980
.340.000000000, 341.000000000, 342.000000000, 343.000000000,
1981
.344.000000000, 345.000000000, 346.000000000, 347.000000000,
1982
.348.000000000, 349.000000000, 350.000000000, 351.000000000,
1983
.353.000000000, 355.000000000, 360.000000000, 365.000000000,
1984
.370.000000000, 375.000000000, 380.000000000, 385.000000000,
1985
.390.000000000, 400.000000000, 410.000000000, 420.000000000,
1986
.430.000000000, 440.000000000, 450.000000000, 460.000000000,
1987
.470.000000000, 480.000000000, 490.000000000, 500.000000000,
1988
.550.000000000, 600.000000000, 650.000000000, 700.000000000,
1989
.750.000000000, 800.000000000, 850.000000000, 900.000000000,
1990
.950.000000000, 1000.000000000, 1100.000000000, 1150.000000000,
1991
.1200.000000000, 1250.000000000, 1300.000000000, 1350.000000000,
1992
.1400.000000000, 1600.000000000, 1800.000000000, 2000.000000000/
1993
*
1994
DATA (yc(i),i=1,120)/
1995
.12.283000000, 12.212000000, 12.118000000, 11.982000000,
1996
.11.780000000, 11.620000000, 11.394000000, 11.333000000,
1997
.11.265000000, 11.191000000, 11.099000000, 10.994000000,
1998
.10.862000000, 10.705000000, 10.510000000, 10.290000000,
1999
.10.112000000, 10.104000000, 10.202000000, 10.238000000,
2000
.10.197000000, 10.115000000, 10.019000000, 9.919000000,
2001
. 9.822000000, 9.726000000, 9.632000000, 9.531000000,
2002
. 9.432000000, 9.334000000, 9.234000000, 9.134000000,
2003
. 9.027000000, 8.907000000, 8.786000000, 8.639000000,
2004
. 8.485000000, 8.351000000, 8.278000000, 8.286000000,
2005
. 8.320000000, 8.336000000, 8.300000000, 8.246000000,
2006
. 8.176000000, 8.107000000, 8.057000000, 7.757000000,
2007
. 7.465000000, 7.044000000, 6.720000000, 6.399000000,
2008
. 6.065000000, 5.815000000, 5.600000000, 5.273000000,
2009
. 5.038000000, 4.795000000, 4.541000000, 4.349000000,
2010
. 4.106000000, 3.978000000, 3.847000000, 3.688000000,
2011
. 3.663000000, 3.589000000, 3.539000000, 3.463000000,
2012
. 3.417000000, 3.317000000, 3.223000000, 3.169000000,
2013
. 3.048000000, 2.449000000, 2.361000000, 2.327000000,
2014
. 2.356000000, 2.470000000, 2.570000000, 2.629000000,
2015
. 2.721000000, 2.834000000, 3.049000000, 3.300000000,
2016
. 3.405000000, 3.492000000, 3.589000000, 3.521000000,
2017
. 3.463000000, 3.260000000, 3.020000000, 2.680000000,
2018
. 2.211000000, 1.710000000, 1.139000000, 0.444000000,
2019
.-0.223000000, -1.060000000, -1.740000000, -2.542000000,
2020
.-7.281000000, -12.708000000, -18.593000000, -24.904000000,
2021
.-31.227000000, -37.790000000, -44.322000000, -50.768000000,
2022
.-57.516000000, -64.187000000, -80.842000000, -92.865000000,
2023
.-131.639000000, -117.769000000, -109.903000000, -122.847000000,
2024
.-101.490000000, -45.186000000, -8.716000000, 17.884000000/
2025
*
2026
n= 120
2027
DO l=1,n
2028
x(l)= xc(l)
2029
y(l)= yc(l)
2030
ENDDO
2031
2032
*.....Set up tridiagonal system.........................................
2033
* b=diagonal, d=offdiagonal, c=right-hand side
2034
*
2035
d(1)= x(2)-x(1)
2036
c(2)= (y(2)-y(1))/d(1)
2037
DO k= 2,n-1
2038
d(k)= x(k+1)-x(k)
2039
b(k)= TWO*(d(k-1)+d(k))
2040
c(k+1)= (y(k+1)-y(k))/d(k)
2041
c(k)= c(k+1)-c(k)
2042
END DO
2043
*
2044
*.....End conditions. third derivatives at x(1) and x(n) obtained
2045
* from divided differences.......................................
2046
*
2047
b(1)= -d(1)
2048
b(n)= -d(n-1)
2049
c(1)= ZERO
2050
c(n)= ZERO
2051
IF (n.gt.3) THEN
2052
c(1)= c(3)/(x(4)-x(2))-c(2)/(x(3)-x(1))
2053
c(n)= c(n-1)/(x(n)-x(n-2))-c(n-2)/(x(n-1)-x(n-3))
2054
c(1)= c(1)*d(1)*d(1)/(x(4)-x(1))
2055
c(n)= -c(n)*d(n-1)*d(n-1)/(x(n)-x(n-3))
2056
END IF
2057
*
2058
DO k=2,n ! forward elimination
2059
t= d(k-1)/b(k-1)
2060
b(k)= b(k)-t*d(k-1)
2061
c(k)= c(k)-t*c(k-1)
2062
END DO
2063
*
2064
c(n)= c(n)/b(n)
2065
*
2066
* back substitution ( makes c the sigma of text)
2067
*
2068
DO k=n-1,1,-1
2069
c(k)= (c(k)-d(k)*c(k+1))/b(k)
2070
END DO
2071
*
2072
*.....Compute polynomial coefficients...................................
2073
*
2074
b(n)= (y(n)-y(n-1))/d(n-1)+d(n-1)*(c(n-1)+c(n)+c(n))
2075
DO k=1,n-1
2076
b(k)= (y(k+1)-y(k))/d(k)-d(k)*(c(k+1)+c(k)+c(k))
2077
d(k)= (c(k+1)-c(k))/d(k)
2078
c(k)= THREE*c(k)
2079
END DO
2080
c(n)= THREE*c(n)
2081
d(n)= d(n-1)
2082
*
2083
RETURN
2084
*
2085
END
2086
*
2087
*------------------------------------------------------------------------
2088
*
2089
SUBROUTINE S4eval3Single2(u,b,c,d,gdim,f,fp,fpp,fppp)
2090
*
2091
* ---------------------------------------------------------------------------
2092
*
2093
* PURPOSE - Evaluate the cubic spline function
2094
* Seval=y(i)+b(i)*(u-x(i))+c(i)*(u-x(i))**2+d(i)*(u-x(i))**3
2095
* where x(i) <= u < x(i+1)
2096
*
2097
* NOTES- if u<x(1), i=1 is used;if u>x(n), i=n is used
2098
*
2099
REAL*8 u
2100
* abscissa at which the spline is to be evaluated
2101
*
2102
INTEGER j,k,n,l,top,gdim
2103
*
2104
REAL*8 xc(120),yc(120)
2105
REAL*8 x(120),y(120)
2106
REAL*8 b(gdim),c(gdim),d(gdim)
2107
* linear,quadratic,cubic coeff
2108
*
2109
REAL*8 f,fp,fpp,fppp
2110
* function, 1st,2nd,3rd deriv
2111
*
2112
INTEGER i
2113
REAL*8 dx
2114
REAL*8 ZERO, TWO, THREE
2115
PARAMETER(ZERO=0.0, TWO=2.0, THREE=3.0)
2116
SAVE i
2117
DATA i/1/
2118
*
2119
* The grid
2120
*
2121
DATA (xc(l),l=1,120)/
2122
.100.000000000, 110.000000000, 120.000000000, 130.000000000,
2123
.140.000000000, 145.000000000, 150.000000000, 151.000000000,
2124
.152.000000000, 153.000000000, 154.000000000, 155.000000000,
2125
.156.000000000, 157.000000000, 158.000000000, 159.000000000,
2126
.160.000000000, 161.000000000, 162.000000000, 163.000000000,
2127
.164.000000000, 165.000000000, 166.000000000, 167.000000000,
2128
.168.000000000, 169.000000000, 170.000000000, 171.000000000,
2129
.172.000000000, 173.000000000, 174.000000000, 175.000000000,
2130
.176.000000000, 177.000000000, 178.000000000, 179.000000000,
2131
.180.000000000, 181.000000000, 182.000000000, 183.000000000,
2132
.184.000000000, 185.000000000, 186.000000000, 187.000000000,
2133
.188.000000000, 189.000000000, 190.000000000, 195.000000000,
2134
.200.000000000, 210.000000000, 220.000000000, 230.000000000,
2135
.240.000000000, 250.000000000, 260.000000000, 270.000000000,
2136
.280.000000000, 290.000000000, 300.000000000, 310.000000000,
2137
.320.000000000, 325.000000000, 330.000000000, 335.000000000,
2138
.336.000000000, 337.000000000, 338.000000000, 339.000000000,
2139
.340.000000000, 341.000000000, 342.000000000, 343.000000000,
2140
.344.000000000, 345.000000000, 346.000000000, 347.000000000,
2141
.348.000000000, 349.000000000, 350.000000000, 351.000000000,
2142
.353.000000000, 355.000000000, 360.000000000, 365.000000000,
2143
.370.000000000, 375.000000000, 380.000000000, 385.000000000,
2144
.390.000000000, 400.000000000, 410.000000000, 420.000000000,
2145
.430.000000000, 440.000000000, 450.000000000, 460.000000000,
2146
.470.000000000, 480.000000000, 490.000000000, 500.000000000,
2147
.550.000000000, 600.000000000, 650.000000000, 700.000000000,
2148
.750.000000000, 800.000000000, 850.000000000, 900.000000000,
2149
.950.000000000, 1000.000000000, 1100.000000000, 1150.000000000,
2150
.1200.000000000, 1250.000000000, 1300.000000000, 1350.000000000,
2151
.1400.000000000, 1600.000000000, 1800.000000000, 2000.000000000/
2152
*
2153
DATA (yc(l),l=1,120)/
2154
.12.283000000, 12.212000000, 12.118000000, 11.982000000,
2155
.11.780000000, 11.620000000, 11.394000000, 11.333000000,
2156
.11.265000000, 11.191000000, 11.099000000, 10.994000000,
2157
.10.862000000, 10.705000000, 10.510000000, 10.290000000,
2158
.10.112000000, 10.104000000, 10.202000000, 10.238000000,
2159
.10.197000000, 10.115000000, 10.019000000, 9.919000000,
2160
. 9.822000000, 9.726000000, 9.632000000, 9.531000000,
2161
. 9.432000000, 9.334000000, 9.234000000, 9.134000000,
2162
. 9.027000000, 8.907000000, 8.786000000, 8.639000000,
2163
. 8.485000000, 8.351000000, 8.278000000, 8.286000000,
2164
. 8.320000000, 8.336000000, 8.300000000, 8.246000000,
2165
. 8.176000000, 8.107000000, 8.057000000, 7.757000000,
2166
. 7.465000000, 7.044000000, 6.720000000, 6.399000000,
2167
. 6.065000000, 5.815000000, 5.600000000, 5.273000000,
2168
. 5.038000000, 4.795000000, 4.541000000, 4.349000000,
2169
. 4.106000000, 3.978000000, 3.847000000, 3.688000000,
2170
. 3.663000000, 3.589000000, 3.539000000, 3.463000000,
2171
. 3.417000000, 3.317000000, 3.223000000, 3.169000000,
2172
. 3.048000000, 2.449000000, 2.361000000, 2.327000000,
2173
. 2.356000000, 2.470000000, 2.570000000, 2.629000000,
2174
. 2.721000000, 2.834000000, 3.049000000, 3.300000000,
2175
. 3.405000000, 3.492000000, 3.589000000, 3.521000000,
2176
. 3.463000000, 3.260000000, 3.020000000, 2.680000000,
2177
. 2.211000000, 1.710000000, 1.139000000, 0.444000000,
2178
.-0.223000000, -1.060000000, -1.740000000, -2.542000000,
2179
.-7.281000000, -12.708000000, -18.593000000, -24.904000000,
2180
.-31.227000000, -37.790000000, -44.322000000, -50.768000000,
2181
.-57.516000000, -64.187000000, -80.842000000, -92.865000000,
2182
.-131.639000000, -117.769000000, -109.903000000, -122.847000000,
2183
.-101.490000000, -45.186000000, -8.716000000, 17.884000000/
2184
*
2185
n= 120
2186
DO l=1,n
2187
x(l)= xc(l)
2188
y(l)= yc(l)
2189
ENDDO
2190
*
2191
*.....First check if u is in the same interval found on the
2192
* last call to Seval.............................................
2193
*
2194
IF ( (i.lt.1) .OR. (i.ge.n) ) i=1
2195
IF ( (u.lt.x(i)) .OR. (u.ge.x(i+1)) ) THEN
2196
i=1
2197
*
2198
* binary search
2199
*
2200
j= n+1
2201
2468 k= (i+j)/2
2202
IF (u.lt.x(k)) THEN
2203
j= k
2204
ELSE
2205
i= k
2206
ENDIF
2207
IF (j.le.i+1) GOTO 2469
2208
GOTO 2468
2209
2469 CONTINUE
2210
ENDIF
2211
*
2212
dx= u-x(i)
2213
*
2214
* evaluate the spline
2215
*
2216
f= y(i)+dx*(b(i)+dx*(c(i)+dx*d(i)))
2217
fp= b(i)+dx*(TWO*c(i) + dx*THREE*d(i))
2218
fpp= TWO*c(i) + dx*SIX*d(i)
2219
fppp= SIX*d(i)
2220
*
2221
RETURN
2222
*
2223
END
2224
Loggerhead is a web-based interface for Breezy
Version: 2.0.1