Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/dcel2.ssp
There are 2 other files named dcel2.ssp in the archive. Click here to see a list.
C                                                                       DCE2  10
C     ..................................................................DCE2  20
C                                                                       DCE2  30
C        SUBROUTINE DCEL2                                               DCE2  40
C                                                                       DCE2  50
C        PURPOSE                                                        DCE2  60
C           COMPUTES THE GENERALIZED COMPLETE ELLIPTIC INTEGRAL OF      DCE2  70
C           SECOND KIND.                                                DCE2  80
C                                                                       DCE2  90
C        USAGE                                                          DCE2 100
C           CALL DCEL2(RES,AK,A,B,IER)                                  DCE2 110
C                                                                       DCE2 120
C        DESCRIPTION OF PARAMETERS                                      DCE2 130
C           RES   - RESULT VALUE IN DOUBLE PRECISION                    DCE2 140
C           AK    - MODULUS (INPUT) IN DOUBLE PRECISION                 DCE2 150
C           A     - DOUBLE PRECISION CONSTANT TERM IN NUMERATOR         DCE2 160
C           B     - DOUBLE PRECISION FACTOR OF QUADRATIC TERM           DCE2 170
C                   IN NUMERATOR                                        DCE2 180
C           IER   - RESULTANT ERROR CODE WHERE                          DCE2 190
C                   IER=0  NO ERROR                                     DCE2 200
C                   IER=1  AK NOT IN RANGE -1 TO +1                     DCE2 210
C                                                                       DCE2 220
C        REMARKS                                                        DCE2 230
C           FOR ABS(AK) GE 1 THE RESULT IS SET TO 1.E75 IF B IS         DCE2 240
C           POSITIVE, TO -1.7D38 IF B IS NEGATIVE.                       DCE2 250
C           SPECIAL CASES ARE                                           DCE2 260
C           K(K) OBTAINED WITH A = 1, B = 1                             DCE2 270
C           E(K) OBTAINED WITH A = 1, B = CK*CK WHERE CK IS             DCE2 280
C           COMPLEMENTARY MODULUS.                                      DCE2 290
C           B(K) OBTAINED WITH A = 1, B = 0                             DCE2 300
C           D(K) OBTAINED WITH A = 0, B = 1                             DCE2 310
C           WHERE K, E, B, D DEFINE SPECIAL CASES OF THE GENERALIZED    DCE2 320
C           COMPLETE ELLIPTIC INTEGRAL OF SECOND KIND IN THE USUAL      DCE2 330
C           NOTATION, AND THE ARGUMENT K OF THESE FUNCTIONS MEANS       DCE2 340
C           THE MODULUS.                                                DCE2 350
C                                                                       DCE2 360
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  DCE2 370
C           NONE                                                        DCE2 380
C                                                                       DCE2 390
C        METHOD                                                         DCE2 400
C           DEFINITION                                                  DCE2 410
C           RES=INTEGRAL((A+B*T*T)/(SQRT((1+T*T)*(1+(CK*T)**2))*(1+T*T))DCE2 420
C           SUMMED OVER T FROM 0 TO INFINITY).                          DCE2 430
C           EVALUATION                                                  DCE2 440
C           LANDENS TRANSFORMATION IS USED FOR CALCULATION.             DCE2 450
C           REFERENCE                                                   DCE2 460
C           R.BULIRSCH, 'NUMERICAL CALCULATION OF ELLIPTIC INTEGRALS    DCE2 470
C           AND ELLIPTIC FUNCTIONS', HANDBOOK SERIES SPECIAL FUNCTIONS, DCE2 480
C           NUMERISCHE MATHEMATIK VOL. 7, 1965, PP. 78-90.              DCE2 490
C                                                                       DCE2 500
C     ..................................................................DCE2 510
C                                                                       DCE2 520
      SUBROUTINE DCEL2(RES,AK,A,B,IER)                                  DCE2 530
      DOUBLE PRECISION RES,AK,A,B,GEO,ARI,AARI,B0,A1                    DCE2 540
      IER=0                                                             DCE2 550
      ARI=2.D0                                                          DCE2 560
      GEO=(0.5D0-AK)+0.5D0                                              DCE2 570
      GEO=GEO+GEO*AK                                                    DCE2 580
      RES=A                                                             DCE2 590
      A1=A+B                                                            DCE2 600
      B0=B+B                                                            DCE2 610
      IF(GEO)1,2,6                                                      DCE2 620
    1 IER=1                                                             DCE2 630
    2 IF(B)3,8,4                                                        DCE2 640
    3 RES=-1.7D38                                                        DCE2 650
      RETURN                                                            DCE2 660
    4 RES=1.7D38                                                         DCE2 670
      RETURN                                                            DCE2 680
    5 GEO=GEO*AARI                                                      DCE2 690
    6 GEO=DSQRT(GEO)                                                    DCE2 700
      GEO=GEO+GEO                                                       DCE2 710
      AARI=ARI                                                          DCE2 720
      ARI=ARI+GEO                                                       DCE2 730
      B0=B0+RES*GEO                                                     DCE2 740
      RES=A1                                                            DCE2 750
      B0=B0+B0                                                          DCE2 760
      A1=B0/ARI+A1                                                      DCE2 770
      IF(GEO/AARI-0.999999995D0)5,7,7                                   DCE2 780
    7 RES=A1/ARI                                                        DCE2 790
      RES=RES+0.57079632679489662D0*RES                                 DCE2 800
    8 RETURN                                                            DCE2 810
      END                                                               DCE2 820