Google
 

Trailing-Edge - PDP-10 Archives - decuslib10-02 - 43,50145/dcel1.ssp
There are 2 other files named dcel1.ssp in the archive. Click here to see a list.
C                                                                       DCE1  10
C     ..................................................................DCE1  20
C                                                                       DCE1  30
C        SUBROUTINE DCEL1                                               DCE1  40
C                                                                       DCE1  50
C        PURPOSE                                                        DCE1  60
C           CALCULATE COMPLETE ELLIPTIC INTEGRAL OF FIRST KIND          DCE1  70
C                                                                       DCE1  80
C        USAGE                                                          DCE1  90
C           CALL DCEL1(RES,AK,IER)                                      DCE1 100
C                                                                       DCE1 110
C        DESCRIPTION OF PARAMETERS                                      DCE1 120
C           RES   - RESULT VALUE IN DOUBLE PRECISION                    DCE1 130
C           AK    - MODULUS (INPUT) IN DOUBLE PRECISION                 DCE1 140
C           IER   - RESULTANT ERROR CODE WHERE                          DCE1 150
C                   IER=0  NO ERROR                                     DCE1 160
C                   IER=1  AK NOT IN RANGE -1 TO +1                     DCE1 170
C                                                                       DCE1 180
C        REMARKS                                                        DCE1 190
C           THE RESULT IS SET TO 1.E75 IF ABS(AK) GE 1                  DCE1 200
C           FOR MODULUS AK AND COMPLEMENTARY MODULUS CK,                DCE1 210
C           EQUATION AK*AK+CK*CK=1.D0 IS USED.                          DCE1 220
C           AK MUST BE IN THE RANGE -1 TO +1                            DCE1 230
C                                                                       DCE1 240
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  DCE1 250
C           NONE                                                        DCE1 260
C                                                                       DCE1 270
C        METHOD                                                         DCE1 280
C           DEFINITION                                                  DCE1 290
C           CEL1(AK)=INTEGRAL(1/SQRT((1+T*T)*(1+(CK*T)**2)), SUMMED     DCE1 300
C           OVER T FROM 0 TO INFINITY).                                 DCE1 310
C           EQUIVALENT ARE THE DEFINITIONS                              DCE1 320
C           CEL1(AK)=INTEGRAL(1/(COS(T)SQRT(1+(CK*TAN(T))**2)),SUMMED   DCE1 330
C           OVER T FROM 0 TO PI/2),                                     DCE1 340
C           CEL1(AK)=INTEGRAL(1/SQRT(1-(AK*SIN(T))**2),SUMMED OVER T    DCE1 350
C           FROM 0 TO PI/2), WHERE K=SQRT(1.-CK*CK).                    DCE1 360
C           EVALUATION                                                  DCE1 370
C           LANDENS TRANSFORMATION IS USED FOR CALCULATION.             DCE1 380
C           REFERENCE                                                   DCE1 390
C           R.BULIRSCH, 'NUMERICAL CALCULATION OF ELLIPTIC INTEGRALS    DCE1 400
C           AND ELLIPTIC FUNCTIONS', HANDBOOK SERIES SPECIAL FUNCTIONS, DCE1 410
C           NUMERISCHE MATHEMATIK VOL. 7, 1965, PP. 78-90.              DCE1 420
C                                                                       DCE1 430
C     ..................................................................DCE1 440
C                                                                       DCE1 450
      SUBROUTINE DCEL1(RES,AK,IER)                                      DCE1 460
      DOUBLE PRECISION RES,AK,GEO,ARI,AARI                              DCE1 470
      IER=0                                                             DCE1 480
      ARI=2.D0                                                          DCE1 490
      GEO=(0.5D0-AK)+0.5D0                                              DCE1 500
      GEO=GEO+GEO*AK                                                    DCE1 510
      RES=0.5D0                                                         DCE1 520
      IF(GEO)1,2,4                                                      DCE1 530
    1 IER=1                                                             DCE1 540
    2 RES=1.7D38                                                         DCE1 550
      RETURN                                                            DCE1 560
    3 GEO=GEO*AARI                                                      DCE1 570
    4 GEO=DSQRT(GEO)                                                    DCE1 580
      GEO=GEO+GEO                                                       DCE1 590
      AARI=ARI                                                          DCE1 600
      ARI=ARI+GEO                                                       DCE1 610
      RES=RES+RES                                                       DCE1 620
      IF(GEO/AARI-0.999999995D0)3,5,5                                   DCE1 630
    5 RES=RES/ARI*6.2831853071795865D0                                  DCE1 640
      RETURN                                                            DCE1 650
      END                                                               DCE1 660