github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/mathext/internal/amos/amoslib/zacon.f (about) 1 SUBROUTINE ZACON(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, FNUL, 2 * TOL, ELIM, ALIM) 3 C***BEGIN PROLOGUE ZACON 4 C***REFER TO ZBESK,ZBESH 5 C 6 C ZACON APPLIES THE ANALYTIC CONTINUATION FORMULA 7 C 8 C K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN) 9 C MP=PI*MR*CMPLX(0.0,1.0) 10 C 11 C TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT 12 C HALF Z PLANE 13 C 14 C***ROUTINES CALLED ZBINU,ZBKNU,ZS1S2,D1MACH,ZABS,ZMLT 15 C***END PROLOGUE ZACON 16 C COMPLEX CK,CONE,CSCL,CSCR,CSGN,CSPN,CY,CZERO,C1,C2,RZ,SC1,SC2,ST, 17 C *S1,S2,Y,Z,ZN 18 DOUBLE PRECISION ALIM, ARG, ASCLE, AS2, AZN, BRY, BSCLE, CKI, 19 * CKR, CONER, CPN, CSCL, CSCR, CSGNI, CSGNR, CSPNI, CSPNR, 20 * CSR, CSRR, CSSR, CYI, CYR, C1I, C1M, C1R, C2I, C2R, ELIM, FMR, 21 * FN, FNU, FNUL, PI, PTI, PTR, RAZN, RL, RZI, RZR, SC1I, SC1R, 22 * SC2I, SC2R, SGN, SPN, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, 23 * YY, ZEROR, ZI, ZNI, ZNR, ZR, D1MACH, ZABS 24 INTEGER I, INU, IUF, KFLAG, KODE, MR, N, NN, NW, NZ 25 DIMENSION YR(N), YI(N), CYR(2), CYI(2), CSSR(3), CSRR(3), BRY(3) 26 DATA PI / 3.14159265358979324D0 / 27 DATA ZEROR,CONER / 0.0D0,1.0D0 / 28 NZ = 0 29 ZNR = -ZR 30 ZNI = -ZI 31 NN = N 32 CALL ZBINU(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, FNUL, TOL, 33 * ELIM, ALIM) 34 IF (NW.LT.0) GO TO 90 35 C----------------------------------------------------------------------- 36 C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION 37 C----------------------------------------------------------------------- 38 NN = MIN0(2,N) 39 CALL ZBKNU(ZNR, ZNI, FNU, KODE, NN, CYR, CYI, NW, TOL, ELIM, ALIM) 40 IF (NW.NE.0) GO TO 90 41 S1R = CYR(1) 42 S1I = CYI(1) 43 FMR = DBLE(FLOAT(MR)) 44 SGN = -DSIGN(PI,FMR) 45 CSGNR = ZEROR 46 CSGNI = SGN 47 IF (KODE.EQ.1) GO TO 10 48 YY = -ZNI 49 CPN = DCOS(YY) 50 SPN = DSIN(YY) 51 CALL ZMLT(CSGNR, CSGNI, CPN, SPN, CSGNR, CSGNI) 52 10 CONTINUE 53 C----------------------------------------------------------------------- 54 C CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE 55 C WHEN FNU IS LARGE 56 C----------------------------------------------------------------------- 57 INU = INT(SNGL(FNU)) 58 ARG = (FNU-DBLE(FLOAT(INU)))*SGN 59 CPN = DCOS(ARG) 60 SPN = DSIN(ARG) 61 CSPNR = CPN 62 CSPNI = SPN 63 IF (MOD(INU,2).EQ.0) GO TO 20 64 CSPNR = -CSPNR 65 CSPNI = -CSPNI 66 20 CONTINUE 67 IUF = 0 68 C1R = S1R 69 C1I = S1I 70 C2R = YR(1) 71 C2I = YI(1) 72 ASCLE = 1.0D+3*D1MACH(1)/TOL 73 IF (KODE.EQ.1) GO TO 30 74 CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF) 75 NZ = NZ + NW 76 SC1R = C1R 77 SC1I = C1I 78 30 CONTINUE 79 CALL ZMLT(CSPNR, CSPNI, C1R, C1I, STR, STI) 80 CALL ZMLT(CSGNR, CSGNI, C2R, C2I, PTR, PTI) 81 YR(1) = STR + PTR 82 YI(1) = STI + PTI 83 IF (N.EQ.1) RETURN 84 CSPNR = -CSPNR 85 CSPNI = -CSPNI 86 S2R = CYR(2) 87 S2I = CYI(2) 88 C1R = S2R 89 C1I = S2I 90 C2R = YR(2) 91 C2I = YI(2) 92 IF (KODE.EQ.1) GO TO 40 93 CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF) 94 NZ = NZ + NW 95 SC2R = C1R 96 SC2I = C1I 97 40 CONTINUE 98 CALL ZMLT(CSPNR, CSPNI, C1R, C1I, STR, STI) 99 CALL ZMLT(CSGNR, CSGNI, C2R, C2I, PTR, PTI) 100 YR(2) = STR + PTR 101 YI(2) = STI + PTI 102 IF (N.EQ.2) RETURN 103 CSPNR = -CSPNR 104 CSPNI = -CSPNI 105 AZN = ZABS(CMPLX(ZNR,ZNI,kind=KIND(1.0D0))) 106 RAZN = 1.0D0/AZN 107 STR = ZNR*RAZN 108 STI = -ZNI*RAZN 109 RZR = (STR+STR)*RAZN 110 RZI = (STI+STI)*RAZN 111 FN = FNU + 1.0D0 112 CKR = FN*RZR 113 CKI = FN*RZI 114 C----------------------------------------------------------------------- 115 C SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS 116 C----------------------------------------------------------------------- 117 CSCL = 1.0D0/TOL 118 CSCR = TOL 119 CSSR(1) = CSCL 120 CSSR(2) = CONER 121 CSSR(3) = CSCR 122 CSRR(1) = CSCR 123 CSRR(2) = CONER 124 CSRR(3) = CSCL 125 BRY(1) = ASCLE 126 BRY(2) = 1.0D0/ASCLE 127 BRY(3) = D1MACH(2) 128 AS2 = ZABS(CMPLX(S2R,S2I,kind=KIND(1.0D0))) 129 KFLAG = 2 130 IF (AS2.GT.BRY(1)) GO TO 50 131 KFLAG = 1 132 GO TO 60 133 50 CONTINUE 134 IF (AS2.LT.BRY(2)) GO TO 60 135 KFLAG = 3 136 60 CONTINUE 137 BSCLE = BRY(KFLAG) 138 S1R = S1R*CSSR(KFLAG) 139 S1I = S1I*CSSR(KFLAG) 140 S2R = S2R*CSSR(KFLAG) 141 S2I = S2I*CSSR(KFLAG) 142 CSR = CSRR(KFLAG) 143 DO 80 I=3,N 144 STR = S2R 145 STI = S2I 146 S2R = CKR*STR - CKI*STI + S1R 147 S2I = CKR*STI + CKI*STR + S1I 148 S1R = STR 149 S1I = STI 150 C1R = S2R*CSR 151 C1I = S2I*CSR 152 STR = C1R 153 STI = C1I 154 C2R = YR(I) 155 C2I = YI(I) 156 IF (KODE.EQ.1) GO TO 70 157 IF (IUF.LT.0) GO TO 70 158 CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF) 159 NZ = NZ + NW 160 SC1R = SC2R 161 SC1I = SC2I 162 SC2R = C1R 163 SC2I = C1I 164 IF (IUF.NE.3) GO TO 70 165 IUF = -4 166 S1R = SC1R*CSSR(KFLAG) 167 S1I = SC1I*CSSR(KFLAG) 168 S2R = SC2R*CSSR(KFLAG) 169 S2I = SC2I*CSSR(KFLAG) 170 STR = SC2R 171 STI = SC2I 172 70 CONTINUE 173 PTR = CSPNR*C1R - CSPNI*C1I 174 PTI = CSPNR*C1I + CSPNI*C1R 175 YR(I) = PTR + CSGNR*C2R - CSGNI*C2I 176 YI(I) = PTI + CSGNR*C2I + CSGNI*C2R 177 CKR = CKR + RZR 178 CKI = CKI + RZI 179 CSPNR = -CSPNR 180 CSPNI = -CSPNI 181 IF (KFLAG.GE.3) GO TO 80 182 PTR = DABS(C1R) 183 PTI = DABS(C1I) 184 C1M = DMAX1(PTR,PTI) 185 IF (C1M.LE.BSCLE) GO TO 80 186 KFLAG = KFLAG + 1 187 BSCLE = BRY(KFLAG) 188 S1R = S1R*CSR 189 S1I = S1I*CSR 190 S2R = STR 191 S2I = STI 192 S1R = S1R*CSSR(KFLAG) 193 S1I = S1I*CSSR(KFLAG) 194 S2R = S2R*CSSR(KFLAG) 195 S2I = S2I*CSSR(KFLAG) 196 CSR = CSRR(KFLAG) 197 80 CONTINUE 198 RETURN 199 90 CONTINUE 200 NZ = -1 201 IF(NW.EQ.(-2)) NZ=-2 202 RETURN 203 END