Last change
on this file since 4068 was
2,
checked in by lmdz, 25 years ago
|
Initial revision
|
-
Property svn:eol-style set to
native
-
Property svn:executable set to
*
-
Property svn:keywords set to
Author Date Id Revision
|
File size:
1.6 KB
|
Line | |
---|
1 | subroutine spline(x,y,n,yp1,ypn,y2) |
---|
2 | |
---|
3 | c |
---|
4 | |
---|
5 | c Routine to set up the interpolating function for a cubic spline |
---|
6 | |
---|
7 | c interpolation (see "Numerical Recipes" for details). |
---|
8 | |
---|
9 | c |
---|
10 | implicit real (a-h,o-z) |
---|
11 | implicit integer (i-n) |
---|
12 | |
---|
13 | parameter(nllm=4096) |
---|
14 | |
---|
15 | dimension x(n),y(n),y2(n),u(nllm) |
---|
16 | |
---|
17 | c |
---|
18 | c write(6,*)(x(i),i=1,n) |
---|
19 | c write(6,*)(y(i),i=1,n) |
---|
20 | |
---|
21 | if(yp1.gt.0.99E30) then |
---|
22 | c the lower boundary condition is set |
---|
23 | y2(1)=0. |
---|
24 | c either to be "natural" |
---|
25 | u(1)=0. |
---|
26 | |
---|
27 | else |
---|
28 | c or else to have a specified first |
---|
29 | y2(1)=-0.5 |
---|
30 | c derivative |
---|
31 | u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) |
---|
32 | |
---|
33 | end if |
---|
34 | |
---|
35 | do 11 i=2,n-1 |
---|
36 | c decomposition loop of the tridiagonal |
---|
37 | sig=(x(i)-x(i-1))/(x(i+1)-x(i-1)) |
---|
38 | c algorithm. Y2 and U are used |
---|
39 | p=sig*y2(i-1)+2. |
---|
40 | c for temporary storage of the decompo- |
---|
41 | y2(i)=(sig-1.)/p |
---|
42 | c sed factors |
---|
43 | u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) |
---|
44 | |
---|
45 | . /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p |
---|
46 | |
---|
47 | 11 continue |
---|
48 | |
---|
49 | if(ypn.gt.0.99E30) then |
---|
50 | c the upper boundary condition is set |
---|
51 | qn=0. |
---|
52 | c either to be "natural" |
---|
53 | un=0. |
---|
54 | |
---|
55 | else |
---|
56 | c or else to have a specified first |
---|
57 | qn=0.5 |
---|
58 | c derivative |
---|
59 | un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) |
---|
60 | |
---|
61 | end if |
---|
62 | |
---|
63 | y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.) |
---|
64 | |
---|
65 | do 12 k=n-1,1,-1 |
---|
66 | c this is the backsubstitution loop of |
---|
67 | y2(k)=y2(k)*y2(k+1)+u(k) |
---|
68 | c the tridiagonal algorithm |
---|
69 | 12 continue |
---|
70 | |
---|
71 | c |
---|
72 | |
---|
73 | return |
---|
74 | |
---|
75 | end |
---|
76 | |
---|
Note: See
TracBrowser
for help on using the repository browser.