home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / polysptp / dtlep.for < prev    next >
Text File  |  1985-11-29  |  3KB  |  98 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DTLEP
  5. C
  6. C        PURPOSE
  7. C           A SERIES EXPANSION IN LEGENDRE POLYNOMIALS WITH INDEPENDENT
  8. C           VARIABLE X IS TRANSFORMED TO A POLYNOMIAL WITH INDEPENDENT
  9. C           VARIABLE Z, WHERE X=A*Z+B
  10. C
  11. C        USAGE
  12. C           CALL DTLEP(A,B,POL,N,C,WORK)
  13. C
  14. C        DESCRIPTION OF PARAMETERS
  15. C           A     - FACTOR OF LINEAR TERM IN GIVEN LINEAR TRANSFORMATION
  16. C                   DOUBLE PRECISION VARIABLE
  17. C           B     - CONSTANT TERM IN GIVEN LINEAR TRANSFORMATION
  18. C                   DOUBLE PRECISION VARIABLE
  19. C           POL   - COEFFICIENT VECTOR OF POLYNOMIAL (RESULTANT VALUE)
  20. C                   COEFFICIENTS ARE ORDERED FROM LOW TO HIGH
  21. C                   DOUBLE PRECISION VECTOR
  22. C           N     - DIMENSION OF COEFFICIENT VECTORS POL AND C
  23. C           C     - GIVEN COEFFICIENT VECTOR OF EXPANSION
  24. C                   COEFFICIENTS ARE ORDERED FROM LOW TO HIGH
  25. C                   POL AND C MAY BE IDENTICALLY LOCATED
  26. C                   DOUBLE PRECISION VECTOR
  27. C           WORK  - WORKING STORAGE OF DIMENSION 2*N
  28. C                   DOUBLE PRECISION ARRAY
  29. C
  30. C        REMARKS
  31. C           COEFFICIENT VECTOR C REMAINS UNCHANGED IF NOT COINCIDING
  32. C           WITH COEFFICIENT VECTOR POL.
  33. C           OPERATION IS BYPASSED IN CASE N LESS THAN 1.
  34. C           THE LINEAR TRANSFORMATION X=A*Z+B OR Z=(1/A)(X-B) TRANSFORMS
  35. C           THE RANGE (-1,+1) IN X TO THE RANGE (ZL,ZR) IN Z, WHERE
  36. C           ZL=-(1+B)/A AND ZR=(1-B)/A.
  37. C           FOR GIVEN ZL, ZR WE HAVE A=2/(ZR-ZL) AND B=-(ZR+ZL)/(ZR-ZL)
  38. C
  39. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  40. C           NONE
  41. C
  42. C        METHOD
  43. C           THE TRANSFORMATION IS BASED ON THE RECURRENCE EQUATION
  44. C           FOR LEGENDRE POLYNOMIALS P(N,X)
  45. C           P(N+1,X)=2*X*P(N,X)-P(N-1,X)-(X*P(N,X)-P(N-1,X))/(N+1),
  46. C           WHERE THE FIRST TERM IN BRACKETS IS THE INDEX,
  47. C           THE SECOND IS THE ARGUMENT.
  48. C           STARTING VALUES ARE P(0,X)=1, P(1,X)=X.
  49. C           THE TRANSFORMATION IS IMPLICITLY DEFINED BY MEANS OF
  50. C           X=A*Z+B TOGETHER WITH
  51. C           SUM(POL(I)*Z**(I-1), SUMMED OVER I FROM 1 TO N)
  52. C           =SUM(C(I)*P(I-1,X), SUMMED OVER I FROM 1 TO N).
  53. C
  54. C     ..................................................................
  55. C
  56.       SUBROUTINE DTLEP(A,B,POL,N,C,WORK)
  57. C
  58.       DIMENSION POL(1),C(1),WORK(1)
  59.       DOUBLE PRECISION A,B,POL,C,WORK,H,P,Q,Q1,FI
  60. C
  61. C        TEST OF DIMENSION
  62.       IF(N-1)2,1,3
  63. C
  64. C        DIMENSION LESS THAN 2
  65.     1 POL(1)=C(1)
  66.     2 RETURN
  67. C
  68.     3 POL(1)=C(1)+B*C(2)
  69.       POL(2)=A*C(2)
  70.       IF(N-2)2,2,4
  71. C
  72. C        INITIALIZATION
  73.     4 WORK(1)=1.D0
  74.       WORK(2)=B
  75.       WORK(3)=0.D0
  76.       WORK(4)=A
  77.       FI=1.D0
  78. C
  79. C        CALCULATE COEFFICIENT VECTOR OF NEXT LEGENDRE POLYNOMIAL
  80. C        AND ADD MULTIPLE OF THIS VECTOR TO POLYNOMIAL POL
  81.       DO 6 J=3,N
  82.       FI=FI+1.D0
  83.       Q=1.D0/FI-1.D0
  84.       Q1=1.D0-Q
  85.       P=0.D0
  86. C
  87.       DO 5 K=2,J
  88.       H=(A*P+B*WORK(2*K-2))*Q1+Q*WORK(2*K-3)
  89.       P=WORK(2*K-2)
  90.       WORK(2*K-2)=H
  91.       WORK(2*K-3)=P
  92.     5 POL(K-1)=POL(K-1)+H*C(J)
  93.       WORK(2*J-1)=0.D0
  94.       WORK(2*J)=A*P*Q1
  95.     6 POL(J)=C(J)*WORK(2*J)
  96.       RETURN
  97.       END
  98.