zacon.f

Go to the documentation of this file.
00001       SUBROUTINE ZACON(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, FNUL,
00002      * TOL, ELIM, ALIM)
00003 C***BEGIN PROLOGUE  ZACON
00004 C***REFER TO  ZBESK,ZBESH
00005 C
00006 C     ZACON APPLIES THE ANALYTIC CONTINUATION FORMULA
00007 C
00008 C         K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
00009 C                 MP=PI*MR*CMPLX(0.0,1.0)
00010 C
00011 C     TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
00012 C     HALF Z PLANE
00013 C
00014 C***ROUTINES CALLED  ZBINU,ZBKNU,ZS1S2,D1MACH,XZABS,ZMLT
00015 C***END PROLOGUE  ZACON
00016 C     COMPLEX CK,CONE,CSCL,CSCR,CSGN,CSPN,CY,CZERO,C1,C2,RZ,SC1,SC2,ST,
00017 C    *S1,S2,Y,Z,ZN
00018       DOUBLE PRECISION ALIM, ARG, ASCLE, AS2, AZN, BRY, BSCLE, CKI,
00019      * CKR, CONER, CPN, CSCL, CSCR, CSGNI, CSGNR, CSPNI, CSPNR,
00020      * CSR, CSRR, CSSR, CYI, CYR, C1I, C1M, C1R, C2I, C2R, ELIM, FMR,
00021      * FN, FNU, FNUL, PI, PTI, PTR, RAZN, RL, RZI, RZR, SC1I, SC1R,
00022      * SC2I, SC2R, SGN, SPN, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR,
00023      * YY, ZEROR, ZI, ZNI, ZNR, ZR, D1MACH, XZABS
00024       INTEGER I, INU, IUF, KFLAG, KODE, MR, N, NN, NW, NZ
00025       DIMENSION YR(N), YI(N), CYR(2), CYI(2), CSSR(3), CSRR(3), BRY(3)
00026       DATA PI / 3.14159265358979324D0 /
00027       DATA ZEROR,CONER / 0.0D0,1.0D0 /
00028       NZ = 0
00029       ZNR = -ZR
00030       ZNI = -ZI
00031       NN = N
00032       CALL ZBINU(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, FNUL, TOL,
00033      * ELIM, ALIM)
00034       IF (NW.LT.0) GO TO 90
00035 C-----------------------------------------------------------------------
00036 C     ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
00037 C-----------------------------------------------------------------------
00038       NN = MIN0(2,N)
00039       CALL ZBKNU(ZNR, ZNI, FNU, KODE, NN, CYR, CYI, NW, TOL, ELIM, ALIM)
00040       IF (NW.NE.0) GO TO 90
00041       S1R = CYR(1)
00042       S1I = CYI(1)
00043       FMR = DBLE(FLOAT(MR))
00044       SGN = -DSIGN(PI,FMR)
00045       CSGNR = ZEROR
00046       CSGNI = SGN
00047       IF (KODE.EQ.1) GO TO 10
00048       YY = -ZNI
00049       CPN = DCOS(YY)
00050       SPN = DSIN(YY)
00051       CALL ZMLT(CSGNR, CSGNI, CPN, SPN, CSGNR, CSGNI)
00052    10 CONTINUE
00053 C-----------------------------------------------------------------------
00054 C     CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
00055 C     WHEN FNU IS LARGE
00056 C-----------------------------------------------------------------------
00057       INU = INT(SNGL(FNU))
00058       ARG = (FNU-DBLE(FLOAT(INU)))*SGN
00059       CPN = DCOS(ARG)
00060       SPN = DSIN(ARG)
00061       CSPNR = CPN
00062       CSPNI = SPN
00063       IF (MOD(INU,2).EQ.0) GO TO 20
00064       CSPNR = -CSPNR
00065       CSPNI = -CSPNI
00066    20 CONTINUE
00067       IUF = 0
00068       C1R = S1R
00069       C1I = S1I
00070       C2R = YR(1)
00071       C2I = YI(1)
00072       ASCLE = 1.0D+3*D1MACH(1)/TOL
00073       IF (KODE.EQ.1) GO TO 30
00074       CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF)
00075       NZ = NZ + NW
00076       SC1R = C1R
00077       SC1I = C1I
00078    30 CONTINUE
00079       CALL ZMLT(CSPNR, CSPNI, C1R, C1I, STR, STI)
00080       CALL ZMLT(CSGNR, CSGNI, C2R, C2I, PTR, PTI)
00081       YR(1) = STR + PTR
00082       YI(1) = STI + PTI
00083       IF (N.EQ.1) RETURN
00084       CSPNR = -CSPNR
00085       CSPNI = -CSPNI
00086       S2R = CYR(2)
00087       S2I = CYI(2)
00088       C1R = S2R
00089       C1I = S2I
00090       C2R = YR(2)
00091       C2I = YI(2)
00092       IF (KODE.EQ.1) GO TO 40
00093       CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF)
00094       NZ = NZ + NW
00095       SC2R = C1R
00096       SC2I = C1I
00097    40 CONTINUE
00098       CALL ZMLT(CSPNR, CSPNI, C1R, C1I, STR, STI)
00099       CALL ZMLT(CSGNR, CSGNI, C2R, C2I, PTR, PTI)
00100       YR(2) = STR + PTR
00101       YI(2) = STI + PTI
00102       IF (N.EQ.2) RETURN
00103       CSPNR = -CSPNR
00104       CSPNI = -CSPNI
00105       AZN = XZABS(ZNR,ZNI)
00106       RAZN = 1.0D0/AZN
00107       STR = ZNR*RAZN
00108       STI = -ZNI*RAZN
00109       RZR = (STR+STR)*RAZN
00110       RZI = (STI+STI)*RAZN
00111       FN = FNU + 1.0D0
00112       CKR = FN*RZR
00113       CKI = FN*RZI
00114 C-----------------------------------------------------------------------
00115 C     SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS
00116 C-----------------------------------------------------------------------
00117       CSCL = 1.0D0/TOL
00118       CSCR = TOL
00119       CSSR(1) = CSCL
00120       CSSR(2) = CONER
00121       CSSR(3) = CSCR
00122       CSRR(1) = CSCR
00123       CSRR(2) = CONER
00124       CSRR(3) = CSCL
00125       BRY(1) = ASCLE
00126       BRY(2) = 1.0D0/ASCLE
00127       BRY(3) = D1MACH(2)
00128       AS2 = XZABS(S2R,S2I)
00129       KFLAG = 2
00130       IF (AS2.GT.BRY(1)) GO TO 50
00131       KFLAG = 1
00132       GO TO 60
00133    50 CONTINUE
00134       IF (AS2.LT.BRY(2)) GO TO 60
00135       KFLAG = 3
00136    60 CONTINUE
00137       BSCLE = BRY(KFLAG)
00138       S1R = S1R*CSSR(KFLAG)
00139       S1I = S1I*CSSR(KFLAG)
00140       S2R = S2R*CSSR(KFLAG)
00141       S2I = S2I*CSSR(KFLAG)
00142       CSR = CSRR(KFLAG)
00143       DO 80 I=3,N
00144         STR = S2R
00145         STI = S2I
00146         S2R = CKR*STR - CKI*STI + S1R
00147         S2I = CKR*STI + CKI*STR + S1I
00148         S1R = STR
00149         S1I = STI
00150         C1R = S2R*CSR
00151         C1I = S2I*CSR
00152         STR = C1R
00153         STI = C1I
00154         C2R = YR(I)
00155         C2I = YI(I)
00156         IF (KODE.EQ.1) GO TO 70
00157         IF (IUF.LT.0) GO TO 70
00158         CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF)
00159         NZ = NZ + NW
00160         SC1R = SC2R
00161         SC1I = SC2I
00162         SC2R = C1R
00163         SC2I = C1I
00164         IF (IUF.NE.3) GO TO 70
00165         IUF = -4
00166         S1R = SC1R*CSSR(KFLAG)
00167         S1I = SC1I*CSSR(KFLAG)
00168         S2R = SC2R*CSSR(KFLAG)
00169         S2I = SC2I*CSSR(KFLAG)
00170         STR = SC2R
00171         STI = SC2I
00172    70   CONTINUE
00173         PTR = CSPNR*C1R - CSPNI*C1I
00174         PTI = CSPNR*C1I + CSPNI*C1R
00175         YR(I) = PTR + CSGNR*C2R - CSGNI*C2I
00176         YI(I) = PTI + CSGNR*C2I + CSGNI*C2R
00177         CKR = CKR + RZR
00178         CKI = CKI + RZI
00179         CSPNR = -CSPNR
00180         CSPNI = -CSPNI
00181         IF (KFLAG.GE.3) GO TO 80
00182         PTR = DABS(C1R)
00183         PTI = DABS(C1I)
00184         C1M = DMAX1(PTR,PTI)
00185         IF (C1M.LE.BSCLE) GO TO 80
00186         KFLAG = KFLAG + 1
00187         BSCLE = BRY(KFLAG)
00188         S1R = S1R*CSR
00189         S1I = S1I*CSR
00190         S2R = STR
00191         S2I = STI
00192         S1R = S1R*CSSR(KFLAG)
00193         S1I = S1I*CSSR(KFLAG)
00194         S2R = S2R*CSSR(KFLAG)
00195         S2I = S2I*CSSR(KFLAG)
00196         CSR = CSRR(KFLAG)
00197    80 CONTINUE
00198       RETURN
00199    90 CONTINUE
00200       NZ = -1
00201       IF(NW.EQ.(-2)) NZ=-2
00202       RETURN
00203       END
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines