QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
qpms/amos/zuoik.f

195 lines
6.4 KiB
Fortran
Raw Permalink Blame History

SUBROUTINE ZUOIK(ZR, ZI, FNU, KODE, IKFLG, N, YR, YI, NUF, TOL,
* ELIM, ALIM)
C***BEGIN PROLOGUE ZUOIK
C***REFER TO ZBESI,ZBESK,ZBESH
C
C ZUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC
C EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM
C (IN LOGARITHMIC FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW
C WHERE ALIM.LT.ELIM. IF THE MAGNITUDE, BASED ON THE LEADING
C EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN
C THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER
C MULTIPLIERS (IN LOGARITHMIC FORM) IS MADE BASED ON ELIM. HERE
C EXP(-ELIM)=SMALLEST MACHINE NUMBER*1.0E+3 AND EXP(-ALIM)=
C EXP(-ELIM)/TOL
C
C IKFLG=1 MEANS THE I SEQUENCE IS TESTED
C =2 MEANS THE K SEQUENCE IS TESTED
C NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE
C =-1 MEANS AN OVERFLOW WOULD OCCUR
C IKFLG=1 AND NUF.GT.0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO
C THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE
C IKFLG=2 AND NUF.EQ.N MEANS ALL Y VALUES WERE SET TO ZERO
C IKFLG=2 AND 0.LT.NUF.LT.N NOT CONSIDERED. Y MUST BE SET BY
C ANOTHER ROUTINE
C
C***ROUTINES CALLED ZUCHK,ZUNHJ,ZUNIK,D1MACH,AZABS,AZLOG
C***END PROLOGUE ZUOIK
C COMPLEX ARG,ASUM,BSUM,CWRK,CZ,CZERO,PHI,SUM,Y,Z,ZB,ZETA1,ZETA2,ZN,
C *ZR
DOUBLE PRECISION AARG, AIC, ALIM, APHI, ARGI, ARGR, ASUMI, ASUMR,
* ASCLE, AX, AY, BSUMI, BSUMR, CWRKI, CWRKR, CZI, CZR, ELIM, FNN,
* FNU, GNN, GNU, PHII, PHIR, RCZ, STR, STI, SUMI, SUMR, TOL, YI,
* YR, ZBI, ZBR, ZEROI, ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZI,
* ZNI, ZNR, ZR, ZRI, ZRR, D1MACH, AZABS
INTEGER I, IDUM, IFORM, IKFLG, INIT, KODE, N, NN, NUF, NW
DIMENSION YR(N), YI(N), CWRKR(16), CWRKI(16)
DATA ZEROR,ZEROI / 0.0D0, 0.0D0 /
DATA AIC / 1.265512123484645396D+00 /
NUF = 0
NN = N
ZRR = ZR
ZRI = ZI
IF (ZR.GE.0.0D0) GO TO 10
ZRR = -ZR
ZRI = -ZI
10 CONTINUE
ZBR = ZRR
ZBI = ZRI
AX = DABS(ZR)*1.7321D0
AY = DABS(ZI)
IFORM = 1
IF (AY.GT.AX) IFORM = 2
GNU = DMAX1(FNU,1.0D0)
IF (IKFLG.EQ.1) GO TO 20
FNN = DBLE(FLOAT(NN))
GNN = FNU + FNN - 1.0D0
GNU = DMAX1(GNN,FNN)
20 CONTINUE
C-----------------------------------------------------------------------
C ONLY THE MAGNITUDE OF ARG AND PHI ARE NEEDED ALONG WITH THE
C REAL PARTS OF ZETA1, ZETA2 AND ZB. NO ATTEMPT IS MADE TO GET
C THE SIGN OF THE IMAGINARY PART CORRECT.
C-----------------------------------------------------------------------
IF (IFORM.EQ.2) GO TO 30
INIT = 0
CALL ZUNIK(ZRR, ZRI, GNU, IKFLG, 1, TOL, INIT, PHIR, PHII,
* ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI)
CZR = -ZETA1R + ZETA2R
CZI = -ZETA1I + ZETA2I
GO TO 50
30 CONTINUE
ZNR = ZRI
ZNI = -ZRR
IF (ZI.GT.0.0D0) GO TO 40
ZNR = -ZNR
40 CONTINUE
CALL ZUNHJ(ZNR, ZNI, GNU, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R,
* ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
CZR = -ZETA1R + ZETA2R
CZI = -ZETA1I + ZETA2I
AARG = AZABS(ARGR,ARGI)
50 CONTINUE
IF (KODE.EQ.1) GO TO 60
CZR = CZR - ZBR
CZI = CZI - ZBI
60 CONTINUE
IF (IKFLG.EQ.1) GO TO 70
CZR = -CZR
CZI = -CZI
70 CONTINUE
APHI = AZABS(PHIR,PHII)
RCZ = CZR
C-----------------------------------------------------------------------
C OVERFLOW TEST
C-----------------------------------------------------------------------
IF (RCZ.GT.ELIM) GO TO 210
IF (RCZ.LT.ALIM) GO TO 80
RCZ = RCZ + DLOG(APHI)
IF (IFORM.EQ.2) RCZ = RCZ - 0.25D0*DLOG(AARG) - AIC
IF (RCZ.GT.ELIM) GO TO 210
GO TO 130
80 CONTINUE
C-----------------------------------------------------------------------
C UNDERFLOW TEST
C-----------------------------------------------------------------------
IF (RCZ.LT.(-ELIM)) GO TO 90
IF (RCZ.GT.(-ALIM)) GO TO 130
RCZ = RCZ + DLOG(APHI)
IF (IFORM.EQ.2) RCZ = RCZ - 0.25D0*DLOG(AARG) - AIC
IF (RCZ.GT.(-ELIM)) GO TO 110
90 CONTINUE
DO 100 I=1,NN
YR(I) = ZEROR
YI(I) = ZEROI
100 CONTINUE
NUF = NN
RETURN
110 CONTINUE
ASCLE = 1.0D+3*D1MACH(1)/TOL
CALL AZLOG(PHIR, PHII, STR, STI, IDUM)
CZR = CZR + STR
CZI = CZI + STI
IF (IFORM.EQ.1) GO TO 120
CALL AZLOG(ARGR, ARGI, STR, STI, IDUM)
CZR = CZR - 0.25D0*STR - AIC
CZI = CZI - 0.25D0*STI
120 CONTINUE
AX = DEXP(RCZ)/TOL
AY = CZI
CZR = AX*DCOS(AY)
CZI = AX*DSIN(AY)
CALL ZUCHK(CZR, CZI, NW, ASCLE, TOL)
IF (NW.NE.0) GO TO 90
130 CONTINUE
IF (IKFLG.EQ.2) RETURN
IF (N.EQ.1) RETURN
C-----------------------------------------------------------------------
C SET UNDERFLOWS ON I SEQUENCE
C-----------------------------------------------------------------------
140 CONTINUE
GNU = FNU + DBLE(FLOAT(NN-1))
IF (IFORM.EQ.2) GO TO 150
INIT = 0
CALL ZUNIK(ZRR, ZRI, GNU, IKFLG, 1, TOL, INIT, PHIR, PHII,
* ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI)
CZR = -ZETA1R + ZETA2R
CZI = -ZETA1I + ZETA2I
GO TO 160
150 CONTINUE
CALL ZUNHJ(ZNR, ZNI, GNU, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R,
* ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
CZR = -ZETA1R + ZETA2R
CZI = -ZETA1I + ZETA2I
AARG = AZABS(ARGR,ARGI)
160 CONTINUE
IF (KODE.EQ.1) GO TO 170
CZR = CZR - ZBR
CZI = CZI - ZBI
170 CONTINUE
APHI = AZABS(PHIR,PHII)
RCZ = CZR
IF (RCZ.LT.(-ELIM)) GO TO 180
IF (RCZ.GT.(-ALIM)) RETURN
RCZ = RCZ + DLOG(APHI)
IF (IFORM.EQ.2) RCZ = RCZ - 0.25D0*DLOG(AARG) - AIC
IF (RCZ.GT.(-ELIM)) GO TO 190
180 CONTINUE
YR(NN) = ZEROR
YI(NN) = ZEROI
NN = NN - 1
NUF = NUF + 1
IF (NN.EQ.0) RETURN
GO TO 140
190 CONTINUE
ASCLE = 1.0D+3*D1MACH(1)/TOL
CALL AZLOG(PHIR, PHII, STR, STI, IDUM)
CZR = CZR + STR
CZI = CZI + STI
IF (IFORM.EQ.1) GO TO 200
CALL AZLOG(ARGR, ARGI, STR, STI, IDUM)
CZR = CZR - 0.25D0*STR - AIC
CZI = CZI - 0.25D0*STI
200 CONTINUE
AX = DEXP(RCZ)/TOL
AY = CZI
CZR = AX*DCOS(AY)
CZI = AX*DSIN(AY)
CALL ZUCHK(CZR, CZI, NW, ASCLE, TOL)
IF (NW.NE.0) GO TO 180
RETURN
210 CONTINUE
NUF = -1
RETURN
END
Reference in New Issue View Git Blame Copy Permalink