GNU Octave: libcruft/amos/cwrsk.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CWRSK(ZR, FNU, KODE, N, Y, NZ, CW, TOL, ELIM, ALIM)
00002 C***BEGIN PROLOGUE  CWRSK
00003 C***REFER TO  CBESI,CBESK
00004 C
00005 C     CWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY
00006 C     NORMALIZING THE I FUNCTION RATIOS FROM CRATI BY THE WRONSKIAN
00007 C
00008 C***ROUTINES CALLED  CBKNU,CRATI,R1MACH
00009 C***END PROLOGUE  CWRSK
00010       COMPLEX CINU, CSCL, CT, CW, C1, C2, RCT, ST, Y, ZR
00011       REAL ACT, ACW, ALIM, ASCLE, ELIM, FNU, S1, S2, TOL, YY
00012       INTEGER I, KODE, N, NW, NZ
00013       DIMENSION Y(N), CW(2)
00014 C-----------------------------------------------------------------------
00015 C     I(FNU+I-1,Z) BY BACKWARD RECURRENCE FOR RATIOS
00016 C     Y(I)=I(FNU+I,Z)/I(FNU+I-1,Z) FROM CRATI NORMALIZED BY THE
00017 C     WRONSKIAN WITH K(FNU,Z) AND K(FNU+1,Z) FROM CBKNU.
00018 C-----------------------------------------------------------------------
00019       NZ = 0
00020       CALL CBKNU(ZR, FNU, KODE, 2, CW, NW, TOL, ELIM, ALIM)
00021       IF (NW.NE.0) GO TO 50
00022       CALL CRATI(ZR, FNU, N, Y, TOL)
00023 C-----------------------------------------------------------------------
00024 C     RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z),
00025 C     R(FNU+J-1,Z)=Y(J),  J=1,...,N
00026 C-----------------------------------------------------------------------
00027       CINU = CMPLX(1.0E0,0.0E0)
00028       IF (KODE.EQ.1) GO TO 10
00029       YY = AIMAG(ZR)
00030       S1 = COS(YY)
00031       S2 = SIN(YY)
00032       CINU = CMPLX(S1,S2)
00033    10 CONTINUE
00034 C-----------------------------------------------------------------------
00035 C     ON LOW EXPONENT MACHINES THE K FUNCTIONS CAN BE CLOSE TO BOTH
00036 C     THE UNDER AND OVERFLOW LIMITS AND THE NORMALIZATION MUST BE
00037 C     SCALED TO PREVENT OVER OR UNDERFLOW. CUOIK HAS DETERMINED THAT
00038 C     THE RESULT IS ON SCALE.
00039 C-----------------------------------------------------------------------
00040       ACW = CABS(CW(2))
00041       ASCLE = 1.0E+3*R1MACH(1)/TOL
00042       CSCL = CMPLX(1.0E0,0.0E0)
00043       IF (ACW.GT.ASCLE) GO TO 20
00044       CSCL = CMPLX(1.0E0/TOL,0.0E0)
00045       GO TO 30
00046    20 CONTINUE
00047       ASCLE = 1.0E0/ASCLE
00048       IF (ACW.LT.ASCLE) GO TO 30
00049       CSCL = CMPLX(TOL,0.0E0)
00050    30 CONTINUE
00051       C1 = CW(1)*CSCL
00052       C2 = CW(2)*CSCL
00053       ST = Y(1)
00054 C-----------------------------------------------------------------------
00055 C     CINU=CINU*(CONJG(CT)/CABS(CT))*(1.0E0/CABS(CT) PREVENTS
00056 C     UNDER- OR OVERFLOW PREMATURELY BY SQUARING CABS(CT)
00057 C-----------------------------------------------------------------------
00058       CT = ZR*(C2+ST*C1)
00059       ACT = CABS(CT)
00060       RCT = CMPLX(1.0E0/ACT,0.0E0)
00061       CT = CONJG(CT)*RCT
00062       CINU = CINU*RCT*CT
00063       Y(1) = CINU*CSCL
00064       IF (N.EQ.1) RETURN
00065       DO 40 I=2,N
00066         CINU = ST*CINU
00067         ST = Y(I)
00068         Y(I) = CINU*CSCL
00069    40 CONTINUE
00070       RETURN
00071    50 CONTINUE
00072       NZ = -1
00073       IF(NW.EQ.(-2)) NZ=-2
00074       RETURN
00075       END