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