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

Go to the documentation of this file.
00001       SUBROUTINE ZWRSK(ZRR, ZRI, FNU, KODE, N, YR, YI, NZ, CWR, CWI,
00002      * TOL, ELIM, ALIM)
00003 C***BEGIN PROLOGUE  ZWRSK
00004 C***REFER TO  ZBESI,ZBESK
00005 C
00006 C     ZWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY
00007 C     NORMALIZING THE I FUNCTION RATIOS FROM ZRATI BY THE WRONSKIAN
00008 C
00009 C***ROUTINES CALLED  D1MACH,ZBKNU,ZRATI,XZABS
00010 C***END PROLOGUE  ZWRSK
00011 C     COMPLEX CINU,CSCL,CT,CW,C1,C2,RCT,ST,Y,ZR
00012       DOUBLE PRECISION ACT, ACW, ALIM, ASCLE, CINUI, CINUR, CSCLR, CTI,
00013      * CTR, CWI, CWR, C1I, C1R, C2I, C2R, ELIM, FNU, PTI, PTR, RACT,
00014      * STI, STR, TOL, YI, YR, ZRI, ZRR, XZABS, D1MACH
00015       INTEGER I, KODE, N, NW, NZ
00016       DIMENSION YR(N), YI(N), CWR(2), CWI(2)
00017 C-----------------------------------------------------------------------
00018 C     I(FNU+I-1,Z) BY BACKWARD RECURRENCE FOR RATIOS
00019 C     Y(I)=I(FNU+I,Z)/I(FNU+I-1,Z) FROM CRATI NORMALIZED BY THE
00020 C     WRONSKIAN WITH K(FNU,Z) AND K(FNU+1,Z) FROM CBKNU.
00021 C-----------------------------------------------------------------------
00022       NZ = 0
00023       CALL ZBKNU(ZRR, ZRI, FNU, KODE, 2, CWR, CWI, NW, TOL, ELIM, ALIM)
00024       IF (NW.NE.0) GO TO 50
00025       CALL ZRATI(ZRR, ZRI, FNU, N, YR, YI, TOL)
00026 C-----------------------------------------------------------------------
00027 C     RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z),
00028 C     R(FNU+J-1,Z)=Y(J),  J=1,...,N
00029 C-----------------------------------------------------------------------
00030       CINUR = 1.0D0
00031       CINUI = 0.0D0
00032       IF (KODE.EQ.1) GO TO 10
00033       CINUR = DCOS(ZRI)
00034       CINUI = DSIN(ZRI)
00035    10 CONTINUE
00036 C-----------------------------------------------------------------------
00037 C     ON LOW EXPONENT MACHINES THE K FUNCTIONS CAN BE CLOSE TO BOTH
00038 C     THE UNDER AND OVERFLOW LIMITS AND THE NORMALIZATION MUST BE
00039 C     SCALED TO PREVENT OVER OR UNDERFLOW. CUOIK HAS DETERMINED THAT
00040 C     THE RESULT IS ON SCALE.
00041 C-----------------------------------------------------------------------
00042       ACW = XZABS(CWR(2),CWI(2))
00043       ASCLE = 1.0D+3*D1MACH(1)/TOL
00044       CSCLR = 1.0D0
00045       IF (ACW.GT.ASCLE) GO TO 20
00046       CSCLR = 1.0D0/TOL
00047       GO TO 30
00048    20 CONTINUE
00049       ASCLE = 1.0D0/ASCLE
00050       IF (ACW.LT.ASCLE) GO TO 30
00051       CSCLR = TOL
00052    30 CONTINUE
00053       C1R = CWR(1)*CSCLR
00054       C1I = CWI(1)*CSCLR
00055       C2R = CWR(2)*CSCLR
00056       C2I = CWI(2)*CSCLR
00057       STR = YR(1)
00058       STI = YI(1)
00059 C-----------------------------------------------------------------------
00060 C     CINU=CINU*(CONJG(CT)/CABS(CT))*(1.0D0/CABS(CT) PREVENTS
00061 C     UNDER- OR OVERFLOW PREMATURELY BY SQUARING CABS(CT)
00062 C-----------------------------------------------------------------------
00063       PTR = STR*C1R - STI*C1I
00064       PTI = STR*C1I + STI*C1R
00065       PTR = PTR + C2R
00066       PTI = PTI + C2I
00067       CTR = ZRR*PTR - ZRI*PTI
00068       CTI = ZRR*PTI + ZRI*PTR
00069       ACT = XZABS(CTR,CTI)
00070       RACT = 1.0D0/ACT
00071       CTR = CTR*RACT
00072       CTI = -CTI*RACT
00073       PTR = CINUR*RACT
00074       PTI = CINUI*RACT
00075       CINUR = PTR*CTR - PTI*CTI
00076       CINUI = PTR*CTI + PTI*CTR
00077       YR(1) = CINUR*CSCLR
00078       YI(1) = CINUI*CSCLR
00079       IF (N.EQ.1) RETURN
00080       DO 40 I=2,N
00081         PTR = STR*CINUR - STI*CINUI
00082         CINUI = STR*CINUI + STI*CINUR
00083         CINUR = PTR
00084         STR = YR(I)
00085         STI = YI(I)
00086         YR(I) = CINUR*CSCLR
00087         YI(I) = CINUI*CSCLR
00088    40 CONTINUE
00089       RETURN
00090    50 CONTINUE
00091       NZ = -1
00092       IF(NW.EQ.(-2)) NZ=-2
00093       RETURN
00094       END