tpc/generator/airy.f

0001
0002       SUBROUTINE AIZ(IFUN,IFAC,X0,Y0,GAIR,GAII,IERRO)
0003 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0004 C COMPUTATION OF THE AIRY FUNCTION AI(Z) OR ITS DERIVATIVE AI'(Z)
0005 C THE CODE USES:
0006 C      1. MACLAURIN SERIES FOR |Y|<3 AND -2.5<X<1.3 (Z=X+I*Y)
0007 C      2. GAUSS-LAGUERRE QUADRATURE  FOR |Z|<15 AND  WHEN
0008 C         MACLAURIN SERIES ARE NOT USED.
0009 C      3. ASYMPTOTIC EXPANSION FOR |Z|>15.
0010 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0011 C  INPUTS:
0012 C    IFUN:
0013 C         * IFUN=1, THE CODE COMPUTES AI(Z)
0014 C         * IFUN=2, THE CODE COMPUTES AI'(Z)
0015 C    IFAC:
0016 C         * IFAC=1, THE CODE COMPUTES  AI(Z) OR AI'(Z)
0017 C         * IFAC=2, THE CODE COMPUTES NORMALIZED AI(Z) OR AI'(Z)
0018 C    X0:   REAL PART OF THE ARGUMENT Z
0019 C    Y0:   IMAGINARY PART OF THE ARGUMENT  Z
0020 C
0021 C  OUTPUTS:
0022 C    GAIR: REAL PART OF AI(Z) OR AI'(Z)
0023 C    GAII: IMAGINARY PART OF AI(Z) OR AI'(Z)
0024 C
0025 C    IERRO: ERROR FLAG
0026 C          * IERRO=0, SUCCESSFUL COMPUTATION
0027 C          * IERRO=1, COMPUTATION OUT OF RANGE
0028 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0029 C          MACHINE DEPENDENT CONSTANTS: FUNCTION D1MACH
0030 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0031 C   ACCURACY:
0032 C
0033 C     1) SCALED AIRY FUNCTIONS:
0034 C        RELATIVE ACCURACY BETTER THAN 10**(-13) EXCEPT CLOSE TO
0035 C        THE ZEROS, WHERE 10**(-13) IS THE ABSOLUTE PRECISION.
0036 C        GRADUAL LOSS OF PRECISION TAKES PLACE FOR |Z|>1000
0037 C        (REACHING 10**(-8) ABSOLUTE ACCURACY FOR |Z| CLOSE
0038 C        TO 10**(6)) IN THE CASE OF PHASE(Z) CLOSE TO PI.
0039 C     2) UNSCALED AIRY FUNCTIONS:
0040 C        THE FUNCTION OVERFLOWS/UNDERFLOWS FOR
0041 C        3/2*|Z|**(3/2)>LOG(OVER).
0042 C        FOR |Z|<30:
0043 C        A) RELATIVE ACCURACY FOR THE MODULUS (EXCEPT AT THE
0044 C           ZEROS) BETTER THAN 10**(-13).
0045 C        B) ABSOLUTE ACCURACY FOR MIN(R(Z),1/R(Z)) BETTER
0046 C           THAN 10**(-13), WHERE R(Z)=REAL(AI)/IMAG(AI)
0047 C           OR R(Z)=REAL(AI')/IMAG(AI').
0048 C        FOR |Z|>30, GRADUAL LOSS OF PRECISION TAKES PLACE
0049 C        AS |Z| INCREASES.
0050 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0051 C     AUTHORS:
0052 C        AMPARO GIL    (U. AUTONOMA DE MADRID, MADRID, SPAIN).
0053 C                      E-MAIL: AMPARO.GIL@UAM.ES
0054 C        JAVIER SEGURA (U. CARLOS III DE MADRID, MADRID, SPAIN).
0055 C                      E-MAIL: JSEGURA@MATH.UC3M.ES
0056 C        NICO M. TEMME (CWI, AMSTERDAM, THE NETHERLANDS).
0057 C                      E-MAIL: NICO.TEMME@CWI.NL
0058 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0059 C    REFERENCES:
0060 C         COMPUTING AIRY FUNCTIONS BY NUMERICAL QUADRATURE.
0061 C         NUMERICAL ALGORITHMS (2001).
0062 C         A. GIL, J. SEGURA, N.M. TEMME
0063 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0064       DOUBLE PRECISION X0,Y0,GAIR,GAII,X,W,XD,WD
0065       DOUBLE PRECISION OVER,UNDER,DL1,DL2,COVER,D1MACH
0066       DOUBLE PRECISION PI,PIHAL,PIH3,PISR,A,ALF,THET,R,TH15,S1,C1,
0067      * R32,FACTO,TH025,S3,C3,F23,PI23,SQRT3,XA,YA,F23R,DF1,DF2,
0068      * S11,C11,DEX,DRE,DIMA,GAR,GAI,C,S,U,V,V0,AR,AI,AR1,AI1,
0069      * RO,COE1,COE2,REX,DFR,DFI,AR11,AI11,PHASE
0070       INTEGER IFUN,IFAC,IERRO,IEXPF,IEXPF2,N
0071       DIMENSION X(25),W(25)
0072       DIMENSION XD(25),WD(25)
0073       COMMON/PARAM1/PI,PIHAL
0074       COMMON/PARAM2/PIH3,PISR,A,ALF
0075       COMMON/PARAM3/THET,R,TH15,S1,C1,R32
0076       COMMON/PARAM4/FACTO,TH025,S3,C3
0077       SAVE X,W
0078       SAVE XD,WD
0079       DATA X,W/.283891417994567679D-1,.170985378860034935D0,
0080      *.435871678341770460D0,.823518257913030858D0,1.33452543254227372D0,
0081      *1.96968293206435071D0,2.72998134002859938D0,3.61662161916100897D0,
0082      *4.63102611052654146D0,5.77485171830547694D0,7.05000568630218682D0,
0083      *8.45866437513237792D0,10.0032955242749393D0,11.6866845947722423D0,
0084      *13.5119659344693551D0,15.4826596959377140D0,17.6027156808069112D0,
0085      *19.8765656022785451D0,22.3091856773962780D0,24.9061720212974207D0,
0086      *27.6738320739497190D0,30.6192963295084111D0,33.7506560850239946D0,
0087      *37.0771349708391198D0,40.6093049694341322D0,.143720408803313866D0,
0088      *.230407559241880881D0,.242253045521327626D0,.203636639103440807D0,
0089      *.143760630622921410D0,.869128834706078120D-1,.4541750018329
0090      * 15883D-1,.206118031206069497D-1,.814278821268606972D-2,.280266
0091      *075663377634D-2,.840337441621719716D-3,.219303732907765020D-3,
0092      *.497401659009257760D-4,.978508095920717661D-5,.166542824603725
0093      *563D-5,.244502736801316287D-6,.308537034236207072D-7,.3332960
0094      *72940112245D-8,.306781892316295828D-9,.239331309885375719D-10,
0095      *.157294707710054952D-11,.864936011664392267D-13,.394819815
0096      *638647111D-14,.148271173082850884D-15,.453390377327054458D-17/
0097       DATA XD,WD/.435079659953445D-1,.205779160144678D0,
0098      *.489916161318751D0,.896390483211727D0,1.42582496737580D0,
0099      *2.07903190767599D0,2.85702335104978D0,3.76102058198275D0,
0100      *4.79246521225895D0,5.95303247470003D0,7.24464710774066D0,
0101      *8.66950223642504D0,10.2300817341775D0,11.9291866622602D0,
0102      *13.7699665302828D0,15.7559563095946D0,17.8911203751898D0,
0103      *20.1799048700978D0,22.6273004064466D0,25.2389175786164D0,
0104      *28.0210785229929D0,30.9809287996116D0,34.1265753192057D0,
0105      *37.4672580871163D0,41.0135664833476D0,.576354557898966D-1,
0106      *.139560003272262D0,.187792315011311D0,.187446935256946D0,
0107      *.150716717316301D0,.101069904453380D0,.575274105486025D-1,
0108      *.280625783448681D-1,.117972164134041D-1,.428701743297432D-2,
0109      *.134857915232883D-2,.367337337105948D-3,.865882267841931D-4,
0110      *.176391622890609D-4,.309929190938078D-5,.468479653648208D-6,
0111      *.607273267228907D-7,.672514812555074D-8,.633469931761606D-9,
0112      *.504938861248542D-10,.338602527895834D-11,.189738532450555D-12,
0113      *.881618802142698D-14,.336676636121976D-15,.104594827170761D-16/
0114 CC CONSTANTS CCCCCCCCCCCCCCCCCCCCCCC
0115       PI=3.1415926535897932385D0
0116       PIHAL=1.5707963267948966192D0
0117       PIH3=4.71238898038469D0
0118       F23=.6666666666666666D0
0119       PI23=2.09439510239320D0
0120       PISR=1.77245385090552D0
0121       SQRT3=1.7320508075688772935D0
0122 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0123       YA=Y0
0124       XA=X0
0125       IERRO=0
0126       IEXPF=0
0127       IEXPF2=0
0128       IF (YA.LT.0.D0) YA=-YA
0129       R=SQRT(XA*XA+YA*YA)
0130       R32=R*SQRT(R)
0131       THET=PHASE(XA,YA)
0132       COVER=2.D0/3.D0*R32*ABS(COS(1.5D0*THET))
0133 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
0134       OVER=D1MACH(2)*1.D-3
0135 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
0136       UNDER=D1MACH(1)*1.D+3
0137       DL1=LOG(OVER)
0138       DL2=-LOG(UNDER)
0139       IF (DL1.GT.DL2) OVER=1/UNDER
0140       IF (IFAC.EQ.1) THEN
0141         IF (COVER.GE.LOG(OVER)) THEN
0142 CCC OVERFLOW/UNDERFLOW PROBLEMS.
0143 CCC   CALCULATION ABORTED
0144           IERRO=1
0145           GAIR=0
0146           GAII=0
0147         ENDIF
0148         IF (COVER.GE.(LOG(OVER)*0.2)) IEXPF2=1
0149       ELSE
0150         IF (COVER.GE.(LOG(OVER)*0.2)) IEXPF=1
0151       ENDIF
0152       IF (IERRO.EQ.0) THEN
0153         IF (IFUN.EQ.1) THEN
0154 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0155 CCC CALCULATION OF AI(Z) CCCCCCCCCCCCCCC
0156 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0157 CCCCCCCCCCCC SERIES, INTEGRALS OR EXPANSIONS CCCCCCCCCCCCCCCCCCCCCCCC
0158           IF ((YA.LT.3.D0).AND.(XA.LT.1.3D0).AND.(XA.GT.-2.5D0)) THEN
0159 CCC SERIES CCC
0160             CALL SERAI(XA,YA,GAR,GAI)
0161             IF (IFAC.EQ.2) THEN
0162               THET=PHASE(XA,YA)
0163               TH15=1.5D0*THET
0164               S1=SIN(TH15)
0165               C1=COS(TH15)
0166               F23R=F23*R32
0167               DF1=F23R*C1
0168               DF2=F23R*S1
0169               S11=SIN(DF2)
0170               C11=COS(DF2)
0171               DEX=EXP(DF1)
0172               DRE=DEX*C11
0173               DIMA=DEX*S11
0174               GAIR=DRE*GAR-DIMA*GAI
0175               GAII=DRE*GAI+DIMA*GAR
0176             ELSE
0177               GAIR=GAR
0178               GAII=GAI
0179               IF (Y0.EQ.0.) GAII=0.D0
0180             ENDIF
0181           ELSE
0182             IF (R.GT.15.D0) THEN
0183 CCC ASYMPTOTIC EXPANSIONS CCC
0184               THET=PHASE(XA,YA)
0185               FACTO=0.5D0/PISR*R**(-0.25D0)
0186               IF (THET.GT.PI23) THEN
0187 CCCCCCCCCCC CONNECTION FORMULAE CCCCCCCCCCCCCCCCCCCCCCCCCC
0188 CCC     N= 1: TRANSFORM Z TO W= U+IV=Z EXP( 2 PI I/3)
0189                 N=1
0190                 C=-0.5D0
0191                 S=N*0.5*SQRT3
0192                 U=XA*C-YA*S
0193                 V=XA*S+YA*C
0194                 V0=V
0195                 IF (V.LT.0.D0) V=-V
0196                 THET=PHASE(U,V)
0197                 TH15=1.5D0*THET
0198                 S1=SIN(TH15)
0199                 C1=COS(TH15)
0200                 TH025=THET*0.25D0
0201                 S3=SIN(TH025)
0202                 C3=COS(TH025)
0203                 CALL EXPAI(AR1,AI1)
0204                 IF (V0.LT.0.D0) AI1=-AI1
0205                 AR=-(C*AR1-S*AI1)
0206                 AI=-(S*AR1+C*AI1)
0207                 IF (IEXPF.EQ.0) THEN
0208                   IF (IEXPF2.EQ.0) THEN
0209 CCC     N=-1: TRANSFORM Z TO W= U+IV=Z EXP(-2 PI I/3)
0210                     N=-1
0211                     C=-0.5D0
0212                     S=N*0.5*SQRT3
0213                     U=XA*C-YA*S
0214                     V=XA*S+YA*C
0215                     V0=V
0216                     IF (V.LT.0.D0) V=-V
0217                     THET=PHASE(U,V)
0218                     TH15=1.5D0*THET
0219                     S1=SIN(TH15)
0220                     C1=COS(TH15)
0221                     TH025=THET*0.25D0
0222                     S3=SIN(TH025)
0223                     C3=COS(TH025)
0224                     CALL EXPAI(AR1,AI1)
0225                     IF (V0.LT.0.D0) AI1=-AI1
0226                     THET=PHASE(XA,YA)
0227                     TH15=1.5D0*THET
0228                     S1=SIN(TH15)
0229                     C1=COS(TH15)
0230                     RO=1.333333333333333D0*R32
0231                     COE1=RO*C1
0232                     COE2=RO*S1
0233                     REX=EXP(COE1)
0234                     DFR=REX*COS(COE2)
0235                     DFI=REX*SIN(COE2)
0236                     AR11=DFR*AR1-DFI*AI1
0237                     AI11=DFR*AI1+DFI*AR1
0238                     GAIR=AR-(C*AR11-S*AI11)
0239                     GAII=AI-(S*AR11+C*AI11)
0240                   ELSE
0241                     THET=PHASE(XA,YA)
0242                     TH15=1.5D0*THET
0243                     S1=SIN(TH15)
0244                     C1=COS(TH15)
0245                     GAIR=AR
0246                     GAII=AI
0247                   ENDIF
0248                 ELSE
0249                   GAIR=AR
0250                   GAII=AI
0251                 ENDIF
0252               ELSE
0253 CCCCCCC  ASYMPTOTIC EXPANSION CCCCCCCCCCCCCCC
0254                 THET=PHASE(XA,YA)
0255                 TH15=1.5D0*THET
0256                 S1=SIN(TH15)
0257                 C1=COS(TH15)
0258                 TH025=THET*0.25D0
0259                 S3=SIN(TH025)
0260                 C3=COS(TH025)
0261                 CALL EXPAI(GAIR,GAII)
0262               ENDIF
0263             ELSE
0264 CCC INTEGRALS
0265               A=0.1666666666666666D0
0266               ALF=-A
0267               FACTO=0.280514117723058D0*R**(-0.25D0)
0268               THET=PHASE(XA,YA)
0269               IF (THET.LE.PIHAL) THEN
0270                 TH15=1.5D0*THET
0271                 S1=SIN(TH15)
0272                 C1=COS(TH15)
0273                 TH025=THET*0.25D0
0274                 S3=SIN(TH025)
0275                 C3=COS(TH025)
0276                 CALL AIRY1(X,W,GAIR,GAII)
0277               ENDIF
0278               IF ((THET.GT.PIHAL).AND.(THET.LE.PI23)) THEN
0279                 TH15=1.5D0*THET
0280                 S1=SIN(TH15)
0281                 C1=COS(TH15)
0282                 TH025=THET*0.25D0
0283                 S3=SIN(TH025)
0284                 C3=COS(TH025)
0285                 CALL AIRY2(X,W,GAIR,GAII)
0286               ENDIF
0287               IF (THET.GT.PI23) THEN
0288                 N=1
0289                 C=-0.5D0
0290                 S=N*0.5*SQRT3
0291                 U=XA*C-YA*S
0292                 V=XA*S+YA*C
0293                 V0=V
0294                 IF (V.LT.0.D0) V=-V
0295                 THET=PHASE(U,V)
0296                 IF (THET.LE.PIHAL) THEN
0297                   TH15=1.5D0*THET
0298                   S1=SIN(TH15)
0299                   C1=COS(TH15)
0300                   TH025=THET*0.25D0
0301                   S3=SIN(TH025)
0302                   C3=COS(TH025)
0303                   CALL AIRY1(X,W,AR1,AI1)
0304                 ENDIF
0305                 IF ((THET.GT.PIHAL).AND.(THET.LE.PI23)) THEN
0306                   TH15=1.5D0*THET
0307                   S1=SIN(TH15)
0308                   C1=COS(TH15)
0309                   TH025=THET*0.25D0
0310                   S3=SIN(TH025)
0311                   C3=COS(TH025)
0312                   CALL AIRY2(X,W,AR1,AI1)
0313                 ENDIF
0314                 IF (V0.LT.0.D0) AI1=-AI1
0315                 AR=-(C*AR1-S*AI1)
0316                 AI=-(S*AR1+C*AI1)
0317                 N=-1
0318                 C=-0.5D0
0319                 S=N*0.5*SQRT3
0320                 U=XA*C-YA*S
0321                 V=XA*S+YA*C
0322                 V0=V
0323                 IF (V.LT.0.D0) V=-V
0324                 THET=PHASE(U,V)
0325                 IF (THET.LE.PIHAL) THEN
0326                   TH15=1.5D0*THET
0327                   S1=SIN(TH15)
0328                   C1=COS(TH15)
0329                   TH025=THET*0.25D0
0330                   S3=SIN(TH025)
0331                   C3=COS(TH025)
0332                   CALL AIRY1(X,W,AR1,AI1)
0333                 ENDIF
0334                 IF ((THET.GT.PIHAL).AND.(THET.LE.PI23)) THEN
0335                   TH15=1.5D0*THET
0336                   S1=SIN(TH15)
0337                   C1=COS(TH15)
0338                   TH025=THET*0.25D0
0339                   S3=SIN(TH025)
0340                   C3=COS(TH025)
0341                   CALL AIRY2(X,W,AR1,AI1)
0342                 ENDIF
0343                 IF (V0.LT.0.D0) AI1=-AI1
0344                 THET=PHASE(XA,YA)
0345                 TH15=1.5D0*THET
0346                 S1=SIN(TH15)
0347                 C1=COS(TH15)
0348                 RO=1.333333333333333D0*R32
0349                 COE1=RO*C1
0350                 COE2=RO*S1
0351                 REX=EXP(COE1)
0352                 DFR=REX*COS(COE2)
0353                 DFI=REX*SIN(COE2)
0354                 AR11=DFR*AR1-DFI*AI1
0355                 AI11=DFR*AI1+DFI*AR1
0356                 GAIR=AR-(C*AR11-S*AI11)
0357                 GAII=AI-(S*AR11+C*AI11)
0358               ENDIF
0359             ENDIF
0360             IF (IFAC.EQ.1) THEN
0361 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0362 CCC CALCULATION OF THE UNSCALED AI(Z) CCCCCCCCC
0363 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0364               F23R=F23*R32
0365               DF1=F23R*C1
0366               DF2=F23R*S1
0367               S11=SIN(DF2)
0368               C11=COS(DF2)
0369               DEX=EXP(-DF1)
0370               DRE=DEX*C11
0371               DIMA=-DEX*S11
0372               GAR=DRE*GAIR-DIMA*GAII
0373               GAI=DRE*GAII+DIMA*GAIR
0374               GAIR=GAR
0375               GAII=GAI
0376               IF (Y0.EQ.0.) GAII=0.D0
0377             ENDIF
0378           ENDIF
0379         ELSE
0380 CCCC CALCULATION OF AI�(Z) CCCCCCCCCCC
0381           ALF=0.1666666666666666D0
0382           FACTO=-0.270898621247918D0*R**0.25D0
0383 CCCCCCCCCCCCCCC SERIES OR INTEGRALS CCCCCCCCCCCCCCCCCCCCCCCCCC
0384           IF ((YA.LT.3.D0).AND.(XA.LT.1.3D0).AND.(XA.GT.-2.5D0)) THEN
0385 CCC SERIES
0386             CALL SERAID(XA,YA,GAR,GAI)
0387             IF (IFAC.EQ.2) THEN
0388               THET=PHASE(XA,YA)
0389               TH15=1.5D0*THET
0390               S1=SIN(TH15)
0391               C1=COS(TH15)
0392               F23R=F23*R32
0393               DF1=F23R*C1
0394               DF2=F23R*S1
0395               S11=SIN(DF2)
0396               C11=COS(DF2)
0397               DEX=EXP(DF1)
0398               DRE=DEX*C11
0399               DIMA=DEX*S11
0400               GAIR=DRE*GAR-DIMA*GAI
0401               GAII=DRE*GAI+DIMA*GAR
0402             ELSE
0403               GAIR=GAR
0404               GAII=GAI
0405               IF (Y0.EQ.0.) GAII=0.D0
0406             ENDIF
0407           ELSE
0408             IF (R.GT.15.D0) THEN
0409 CCC  ASYMPTOTIC EXPANSIONS CCCCCCCCCCCCC
0410               THET=PHASE(XA,YA)
0411               FACTO=0.5D0/PISR*R**0.25D0
0412               IF (THET.GT.PI23) THEN
0413 CCCCCCC CONNECTION FORMULAE CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0414 CCC     N= 1: TRANSFORM Z TO W= U+IV=Z EXP( 2 PI I/3)
0415                 N=1
0416                 C=-0.5D0
0417                 S=N*0.5*SQRT3
0418                 U=XA*C-YA*S
0419                 V=XA*S+YA*C
0420                 V0=V
0421                 IF (V.LT.0.D0) V=-V
0422                 THET=PHASE(U,V)
0423                 TH15=1.5D0*THET
0424                 S1=SIN(TH15)
0425                 C1=COS(TH15)
0426                 TH025=THET*0.25D0
0427                 S3=SIN(TH025)
0428                 C3=COS(TH025)
0429                 CALL EXPAID(AR1,AI1)
0430                 IF (V0.LT.0.D0) AI1=-AI1
0431                 AR=-(C*AR1+S*AI1)
0432                 AI=-(-S*AR1+C*AI1)
0433                 IF (IEXPF.EQ.0) THEN
0434                   IF (IEXPF2.EQ.0) THEN
0435 CCC     N=-1: TRANSFORM Z TO W= U+IV=Z EXP(-2 PI I/3)
0436                     N=-1
0437                     C=-0.5D0
0438                     S=N*0.5*SQRT3
0439                     U=XA*C-YA*S
0440                     V=XA*S+YA*C
0441                     V0=V
0442                     IF (V.LT.0.D0) V=-V
0443                     THET=PHASE(U,V)
0444                     TH15=1.5D0*THET
0445                     S1=SIN(TH15)
0446                     C1=COS(TH15)
0447                     TH025=THET*0.25D0
0448                     S3=SIN(TH025)
0449                     C3=COS(TH025)
0450                     CALL EXPAID(AR1,AI1)
0451                     IF (V0.LT.0.D0) AI1=-AI1
0452                     THET=PHASE(XA,YA)
0453                     TH15=1.5D0*THET
0454                     S1=SIN(TH15)
0455                     C1=COS(TH15)
0456                     RO=1.333333333333333D0*R32
0457                     COE1=RO*C1
0458                     COE2=RO*S1
0459                     REX=EXP(COE1)
0460                     DFR=REX*COS(COE2)
0461                     DFI=REX*SIN(COE2)
0462                     AR11=DFR*AR1-DFI*AI1
0463                     AI11=DFR*AI1+DFI*AR1
0464                     GAIR=AR-(C*AR11+S*AI11)
0465                     GAII=AI-(-S*AR11+C*AI11)
0466                   ELSE
0467                     THET=PHASE(XA,YA)
0468                     TH15=1.5D0*THET
0469                     S1=SIN(TH15)
0470                     C1=COS(TH15)
0471                     GAIR=AR
0472                     GAII=AI
0473                   ENDIF
0474                 ELSE
0475                   GAIR=AR
0476                   GAII=AI
0477                 ENDIF
0478               ELSE
0479                 TH15=1.5D0*THET
0480                 S1=SIN(TH15)
0481                 C1=COS(TH15)
0482                 TH025=THET*0.25D0
0483                 S3=SIN(TH025)
0484                 C3=COS(TH025)
0485                 CALL EXPAID(GAIR,GAII)
0486               ENDIF
0487             ELSE
0488 CCC INTEGRALS CCCCCCCCCCCCCCCC
0489               THET=PHASE(XA,YA)
0490               IF (THET.LE.PIHAL) THEN
0491                 TH15=1.5D0*THET
0492                 S1=SIN(TH15)
0493                 C1=COS(TH15)
0494                 TH025=THET*0.25D0
0495                 S3=SIN(TH025)
0496                 C3=COS(TH025)
0497                 CALL AIRY1D(XD,WD,GAIR,GAII)
0498               ENDIF
0499               IF ((THET.GT.PIHAL).AND.(THET.LE.PI23)) THEN
0500                 TH15=1.5D0*THET
0501                 S1=SIN(TH15)
0502                 C1=COS(TH15)
0503                 TH025=THET*0.25D0
0504                 S3=SIN(TH025)
0505                 C3=COS(TH025)
0506                 CALL AIRY2D(XD,WD,GAIR,GAII)
0507               ENDIF
0508               IF (THET.GT.PI23) THEN
0509                 N=1
0510                 C=-0.5D0
0511                 S=N*0.5*SQRT3
0512                 U=XA*C-YA*S
0513                 V=XA*S+YA*C
0514                 V0=V
0515                 IF (V.LT.0.D0) V=-V
0516                 THET=PHASE(U,V)
0517                 IF (THET.LE.PIHAL) THEN
0518                   TH15=1.5D0*THET
0519                   S1=SIN(TH15)
0520                   C1=COS(TH15)
0521                   TH025=THET*0.25D0
0522                   S3=SIN(TH025)
0523                   C3=COS(TH025)
0524                   CALL AIRY1D(XD,WD,AR1,AI1)
0525                 ENDIF
0526                 IF ((THET.GT.PIHAL).AND.(THET.LE.PI23)) THEN
0527                   TH15=1.5D0*THET
0528                   S1=SIN(TH15)
0529                   C1=COS(TH15)
0530                   TH025=THET*0.25D0
0531                   S3=SIN(TH025)
0532                   C3=COS(TH025)
0533                   CALL AIRY2D(XD,WD,AR1,AI1)
0534                 ENDIF
0535                 IF (V0.LT.0.D0) AI1=-AI1
0536                   AR=-(C*AR1+S*AI1)
0537                   AI=-(-S*AR1+C*AI1)
0538                   N=-1
0539                   C=-0.5D0
0540                   S=N*0.5*SQRT3
0541                   U=XA*C-YA*S
0542                   V=XA*S+YA*C
0543                   V0=V
0544                   IF (V.LT.0.D0) V=-V
0545                   THET=PHASE(U,V)
0546                   IF (THET.LE.PIHAL) THEN
0547                     TH15=1.5D0*THET
0548                     S1=SIN(TH15)
0549                     C1=COS(TH15)
0550                     TH025=THET*0.25D0
0551                     S3=SIN(TH025)
0552                     C3=COS(TH025)
0553                     CALL AIRY1D(XD,WD,AR1,AI1)
0554                   ENDIF
0555                   IF ((THET.GT.PIHAL).AND.(THET.LE.PI23)) THEN
0556                     TH15=1.5D0*THET
0557                     S1=SIN(TH15)
0558                     C1=COS(TH15)
0559                     TH025=THET*0.25D0
0560                     S3=SIN(TH025)
0561                     C3=COS(TH025)
0562                     CALL AIRY2D(XD,WD,AR1,AI1)
0563                   ENDIF
0564                   IF (V0.LT.0.D0) AI1=-AI1
0565                   THET=PHASE(XA,YA)
0566                   TH15=1.5D0*THET
0567                   S1=SIN(TH15)
0568                   C1=COS(TH15)
0569                   RO=1.333333333333333D0*R32
0570                   COE1=RO*C1
0571                   COE2=RO*S1
0572                   REX=EXP(COE1)
0573                   DFR=REX*COS(COE2)
0574                   DFI=REX*SIN(COE2)
0575                   AR11=DFR*AR1-DFI*AI1
0576                   AI11=DFR*AI1+DFI*AR1
0577                   GAIR=AR-(C*AR11+S*AI11)
0578                   GAII=AI-(-S*AR11+C*AI11)
0579                 ENDIF
0580               ENDIF
0581               IF (IFAC.EQ.1) THEN
0582 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0583 CCC CALCULATION OF THE UNSCALED AI'(z) CCCCCCC
0584 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0585                 F23R=F23*R32
0586                 DF1=F23R*C1
0587                 DF2=F23R*S1
0588                 S11=SIN(DF2)
0589                 C11=COS(DF2)
0590                 DEX=EXP(-DF1)
0591                 DRE=DEX*C11
0592                 DIMA=-DEX*S11
0593                 GAR=DRE*GAIR-DIMA*GAII
0594                 GAI=DRE*GAII+DIMA*GAIR
0595                 GAIR=GAR
0596                 GAII=GAI
0597                 IF (Y0.EQ.0) GAII=0.D0
0598               ENDIF
0599             ENDIF
0600           ENDIF
0601         ENDIF
0602         IF (Y0.LT.0.D0) GAII=-GAII
0603         RETURN
0604         END
0605         SUBROUTINE  AIRY1(X,W,GAIR,GAII)
0606 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0607 CCC COMPUTES AI(Z) BY GAUSS-LAGUERRE QUADRATURE IN THE SECTOR C
0608 CCC              0 <= PHASE(Z) <= PI/2                        C
0609 CCC                                                           C
0610 CCC INPUTS:                                                   C
0611 CCC      X,W,      NODES AND WEIGHTS FOR THE GAUSSIAN         C
0612 CCC                QUADRATURE                                 C
0613 CCC OUTPUTS:                                                  C
0614 CCC      GAIR, GAII,  REAL AND IMAGINARY PARTS OF AI(Z)       C
0615 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0616         DOUBLE PRECISION X,W,GAIR,GAII
0617         DOUBLE PRECISION PIH3,PISR,A,ALF,THET,R,TH15,S1,C1,
0618      *  R32,FACTO,TH025,S3,C3,SUMAR,SUMAI,DF1,DF1C1,PHI,PHI6,
0619      *  S2,C2,DMODU,DMODU2,FUNR,FUNI,FAC1,FAC2,PHASE
0620         INTEGER I
0621         DIMENSION X(25),W(25)
0622         COMMON/PARAM2/PIH3,PISR,A,ALF
0623         COMMON/PARAM3/THET,R,TH15,S1,C1,R32
0624         COMMON/PARAM4/FACTO,TH025,S3,C3
0625         SUMAR=0.D0
0626         SUMAI=0.D0
0627         DO 1 I=1,25
0628           DF1=1.5D0*X(I)/R32
0629           DF1C1=DF1*C1
0630           PHI=PHASE(2.D0+DF1C1,DF1*S1)
0631           PHI6=PHI/6.D0
0632           S2=SIN(PHI6)
0633           C2=COS(PHI6)
0634           DMODU=SQRT(4.D0+DF1*DF1+4.D0*DF1C1)
0635           DMODU2=DMODU**ALF
0636           FUNR=DMODU2*C2
0637           FUNI=DMODU2*S2
0638           SUMAR=SUMAR+W(I)*FUNR
0639           SUMAI=SUMAI+W(I)*FUNI
0640  1      CONTINUE
0641         FAC1=FACTO*C3
0642         FAC2=FACTO*S3
0643         GAIR=FAC1*SUMAR+FAC2*SUMAI
0644         GAII=FAC1*SUMAI-FAC2*SUMAR
0645         RETURN
0646         END
0647         SUBROUTINE  AIRY2(X,W,GAIR,GAII)
0648 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0649 CCC COMPUTES AI(Z) BY GAUSS-LAGUERRE QUADRATURE IN THE SECTOR C
0650 CCC              PI/2 < PHASE(Z) <= 2PI/3                     C
0651 CCC                                                           C
0652 CCC INPUTS:                                                   C
0653 CCC      X,W,        NODES AND WEIGHTS FOR THE GAUSSIAN       C
0654 CCC                  QUADRATURE                               C
0655 CCC OUTPUTS:                                                  C
0656 CCC      GAIR, GAII, REAL AND IMAGINARY PARTS OF AI(Z)        C
0657 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0658         DOUBLE PRECISION X,W,GAIR,GAII
0659         DOUBLE PRECISION PIH3,PISR,A,ALF,THET,R,TH15,S1,C1,
0660      *  R32,FACTO,TH025,S3,C3,SUMAR,SUMAI,DF1,DF1C1,PHI,PHI6,
0661      *  S2,C2,DMODU,DMODU2,FUNR,FUNI,FAC1,FAC2,PHASE
0662         DOUBLE PRECISION SQR2,SQR2I,TAU,TGTAU,B,ANG,CTAU,CFAC,CT,ST,
0663      *  SUMR,SUMI,TTAU,BETA
0664         INTEGER I
0665         DIMENSION X(25),W(25)
0666         COMMON/PARAM2/PIH3,PISR,A,ALF
0667         COMMON/PARAM3/THET,R,TH15,S1,C1,R32
0668         COMMON/PARAM4/FACTO,TH025,S3,C3
0669         SQR2=1.41421356237310D0
0670         SQR2I=0.707106781186548D0
0671         TAU=TH15-PIH3*0.5D0
0672         TGTAU=DTAN(TAU)
0673         B=5.D0*A
0674         ANG=TAU*B
0675         CTAU=COS(TAU)
0676         CFAC=CTAU**(-B)
0677         CT=COS(ANG)
0678         ST=SIN(ANG)
0679         SUMR=0.D0
0680         SUMI=0.D0
0681         DO 2 I=1,25
0682           DF1=3.D0*X(I)/(CTAU*R32)
0683           DF1C1=DF1*SQR2I*0.5D0
0684           PHI=PHASE(2.D0-DF1C1,DF1C1)
0685           PHI6=PHI/6.D0
0686           TTAU=X(I)*TGTAU
0687           BETA=PHI6-TTAU
0688           S2=SIN(BETA)
0689           C2=COS(BETA)
0690           DMODU=SQRT(4.D0+DF1*DF1*0.25D0-SQR2*DF1)
0691           DMODU2=DMODU**ALF
0692           FUNR=DMODU2*C2
0693           FUNI=DMODU2*S2
0694           SUMR=SUMR+W(I)*FUNR
0695           SUMI=SUMI+W(I)*FUNI
0696  2      CONTINUE
0697         SUMAR=CFAC*(CT*SUMR-ST*SUMI)
0698         SUMAI=CFAC*(CT*SUMI+ST*SUMR)
0699         FAC1=FACTO*C3
0700         FAC2=FACTO*S3
0701         GAIR=FAC1*SUMAR+FAC2*SUMAI
0702         GAII=FAC1*SUMAI-FAC2*SUMAR
0703         RETURN
0704         END
0705         SUBROUTINE AIRY1D(X,W,GAIR,GAII)
0706 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0707 CCC COMPUTES AI'(Z) BY GAUSS-LAGUERRE QUADRATURE IN THE SECTOR C
0708 CCC              0 <= PHASE(Z) <= PI/2                         C
0709 CCC                                                            C
0710 CCC INPUTS:                                                    C
0711 CCC       X,W,      NODES AND WEIGHTS FOR THE GAUSSIAN         C
0712 CCC                 QUADRATURE                                 C
0713 CCC OUTPUTS:                                                   C
0714 CCC       GAIR,GAII, REAL AND IMAGINARY PARTS OF AI'(Z)        C
0715 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0716         DOUBLE PRECISION X,W,GAIR,GAII
0717         DOUBLE PRECISION PIH3,PISR,A,ALF,THET,R,TH15,S1,C1,
0718      *  R32,FACTO,TH025,S3,C3,SUMAR,SUMAI,DF1,DF1C1,PHI,PHI6,
0719      *  S2,C2,DMODU,DMODU2,FUNR,FUNI,FAC1,FAC2,PHASE
0720         INTEGER I
0721         DIMENSION X(25),W(25)
0722         COMMON/PARAM2/PIH3,PISR,A,ALF
0723         COMMON/PARAM3/THET,R,TH15,S1,C1,R32
0724         COMMON/PARAM4/FACTO,TH025,S3,C3
0725         SUMAR=0.D0
0726         SUMAI=0.D0
0727         DO 3 I=1,25
0728           DF1=1.5D0*X(I)/R32
0729           DF1C1=DF1*C1
0730           PHI=PHASE(2.D0+DF1C1,DF1*S1)
0731           PHI6=-PHI*ALF
0732           S2=SIN(PHI6)
0733           C2=COS(PHI6)
0734           DMODU=SQRT(4.D0+DF1*DF1+4.D0*DF1C1)
0735           DMODU2=DMODU**ALF
0736           FUNR=DMODU2*C2
0737           FUNI=DMODU2*S2
0738           SUMAR=SUMAR+W(I)*FUNR
0739           SUMAI=SUMAI+W(I)*FUNI
0740  3      CONTINUE
0741         FAC1=FACTO*C3
0742         FAC2=FACTO*S3
0743         GAIR=FAC1*SUMAR-FAC2*SUMAI
0744         GAII=FAC1*SUMAI+FAC2*SUMAR
0745         RETURN
0746         END
0747         SUBROUTINE  AIRY2D(X,W,GAIR,GAII)
0748 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0749 CCC COMPUTES AI'(Z) BY GAUSS-LAGUERRE QUADRATURE IN THE SECTOR C
0750 CCC              PI/2 < PHASE(Z) <= 3PI/2                      C
0751 CCC                                                            C
0752 CCC INPUTS:                                                    C
0753 CCC      X,W,   NODES AND WEIGHTS FOR THE GAUSSIAN             C
0754 CCC                   QUADRATURE                               C
0755 CCC OUTPUTS:                                                   C
0756 CCC      GAIR,GAII, REAL AND IMAGINARY PARTS OF AI'(Z)         C
0757 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0758         DOUBLE PRECISION X,W,GAIR,GAII
0759         DOUBLE PRECISION PIH3,PISR,A,ALF,THET,R,TH15,S1,C1,
0760      *  R32,FACTO,TH025,S3,C3,SUMAR,SUMAI,DF1,DF1C1,PHI,PHI6,
0761      *  S2,C2,DMODU,DMODU2,FUNR,FUNI,FAC1,FAC2,PHASE
0762         DOUBLE PRECISION SQR2,SQR2I,TAU,TGTAU,B,ANG,CTAU,CFAC,CT,ST,
0763      *  SUMR,SUMI,TTAU,BETA
0764         INTEGER I
0765         DIMENSION X(25),W(25)
0766         COMMON/PARAM2/PIH3,PISR,A,ALF
0767         COMMON/PARAM3/THET,R,TH15,S1,C1,R32
0768         COMMON/PARAM4/FACTO,TH025,S3,C3
0769         SQR2=1.41421356237310D0
0770         SQR2I=0.707106781186548D0
0771         TAU=TH15-PIH3*0.5D0
0772         TGTAU=DTAN(TAU)
0773         B=7.D0*ALF
0774         ANG=TAU*B
0775         CTAU=COS(TAU)
0776         CFAC=CTAU**(-B)
0777         CT=COS(ANG)
0778         ST=SIN(ANG)
0779         SUMR=0.D0
0780         SUMI=0.D0
0781         DO 4 I=1,25
0782           DF1=3.D0*X(I)/(CTAU*R32)
0783           DF1C1=DF1*SQR2I*0.5D0
0784           PHI=PHASE(2.D0-DF1C1,DF1C1)
0785           PHI6=-PHI/6.D0
0786           TTAU=X(I)*TGTAU
0787           BETA=PHI6-TTAU
0788           S2=SIN(BETA)
0789           C2=COS(BETA)
0790           DMODU=SQRT(4.D0+DF1*DF1*0.25D0-SQR2*DF1)
0791           DMODU2=DMODU**ALF
0792           FUNR=DMODU2*C2
0793           FUNI=DMODU2*S2
0794           SUMR=SUMR+W(I)*FUNR
0795           SUMI=SUMI+W(I)*FUNI
0796  4      CONTINUE
0797         SUMAR=CFAC*(CT*SUMR-ST*SUMI)
0798         SUMAI=CFAC*(CT*SUMI+ST*SUMR)
0799         FAC1=FACTO*C3
0800         FAC2=FACTO*S3
0801         GAIR=FAC1*SUMAR-FAC2*SUMAI
0802         GAII=FAC1*SUMAI+FAC2*SUMAR
0803         RETURN
0804         END
0805         DOUBLE PRECISION FUNCTION PHASE(X,Y)
0806         DOUBLE PRECISION PI,PIHAL,X,Y,AY,P
0807 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0808 CCC  COMPUTES THE PHASE OF Z = X + IY, IN (-PI,PI]
0809 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0810         COMMON/PARAM1/PI,PIHAL
0811         IF ((X.EQ.0).AND.(Y.EQ.0)) THEN
0812           P=0.D0
0813         ELSE
0814           AY=ABS(Y)
0815           IF (X.GE.AY) THEN
0816             P=ATAN(AY/X)
0817           ELSEIF ((X+AY).GE.0.D0) THEN
0818             P=PIHAL-ATAN(X/AY)
0819           ELSE
0820             P=PI+ATAN(AY/X)
0821           ENDIF
0822           IF (Y.LT.0.D0) P=-P
0823         ENDIF
0824         PHASE=P
0825         END
0826         SUBROUTINE FGP(X,Y,EPS,FR,FI,GR,GI)
0827 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0828 C    COMPUTES THE FUNCTIONS F AND G FOR THE SERIES  C
0829 C    OF AI'(Z).                                     C
0830 C    THIS ROUTINE IS CALLED BY SERAID.              C
0831 C                                                   C
0832 C    INPUTS:                                        C
0833 C         X,Y,  REAL AND IMAGINARY PARTS OF Z       C
0834 C         EPS,  PRECISION FOR THE COMPUTATION OF    C
0835 C               THE SERIES                          C
0836 C    OUTPUTS:                                       C
0837 C         FR,FI, REAL AND IMAGINARY PARTS OF F      C
0838 C         GR,GI, REAL AND IMAGINARY PARTS OF G      C
0839 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0840         DOUBLE PRECISION X,Y,EPS,FR,FI,GR,GI
0841         INTEGER A,B,K3
0842         DOUBLE PRECISION X2,Y2,U,V,P,Q,CR,CI,DR,DI
0843         X2=X*X
0844         Y2=Y*Y
0845         K3=0
0846         U=X*(X2-3*Y2)
0847         V=Y*(3*X2-Y2)
0848         CR=0.5D0
0849         CI=0.D0
0850         DR=1.D0
0851         DI=0.D0
0852         FR=0.5D0
0853         FI=0.D0
0854         GR=1.D0
0855         GI=0.D0
0856  70     A=(K3+5)*(K3+3)
0857         B=(K3+1)*(K3+3)
0858         P=(U*CR-V*CI)/A
0859         Q=(V*CR+U*CI)/A
0860         CR=P
0861         CI=Q
0862         P=(U*DR-V*DI)/B
0863         Q=(V*DR+U*DI)/B
0864         DR=P
0865         DI=Q
0866         FR=FR+CR
0867         FI=FI+CI
0868         GR=GR+DR
0869         GI=GI+DI
0870         K3=K3+3
0871         IF ((ABS(CR)+ABS(DR)+ABS(CI)+ABS(DI)).GE.EPS) GOTO 70
0872         U=X2-Y2
0873         V=2.D0*X*Y
0874         P=U*FR-V*FI
0875         Q=U*FI+V*FR
0876         FR=P
0877         FI=Q
0878         RETURN
0879         END
0880         SUBROUTINE FG(X,Y,EPS,FR,FI,GR,GI)
0881 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0882 C    COMPUTES THE FUNCTIONS F AND G IN EXPRESSION   C
0883 C    10.4.2 OF ABRAMOWITZ & STEGUN FOR THE SERIES   C
0884 C    OF AI(Z).                                      C
0885 C    THIS ROUTINE IS CALLED BY SERAI.               C
0886 C                                                   C
0887 C    INPUTS:                                        C
0888 C          X,Y,  REAL AND IMAGINARY PARTS OF Z      C
0889 C          EPS,  PRECISION FOR THE COMPUTATION      C
0890 C                OF THE SERIES.                     C
0891 C    OUTPUTS:                                       C
0892 C          FR,FI, REAL AND IMAGINARY PARTS OF F     C
0893 C          GR,GI, REAL AND IMAGINARY PARTS OF G     C
0894 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0895         INTEGER A,B,K3
0896         DOUBLE PRECISION X2,Y2,U,V,P,Q,CR,CI,DR,DI
0897         DOUBLE PRECISION X,Y,EPS,FR,FI,GR,GI
0898         X2=X*X
0899         Y2=Y*Y
0900         K3=0
0901         U=X*(X2-3.D0*Y2)
0902         V=Y*(3.D0*X2-Y2)
0903         CR=1.D0
0904         CI=0.D0
0905         DR=1.D0
0906         DI=0.D0
0907         FR=1.D0
0908         FI=0.D0
0909         GR=1.D0
0910         GI=0.D0
0911  71     A=(K3+2)*(K3+3)
0912         B=(K3+4)*(K3+3)
0913         P=(U*CR-V*CI)/A
0914         Q=(V*CR+U*CI)/A
0915         CR=P
0916         CI=Q
0917         P=(U*DR-V*DI)/B
0918         Q=(V*DR+U*DI)/B
0919         DR=P
0920         DI=Q
0921         FR=FR+CR
0922         FI=FI+CI
0923         GR=GR+DR
0924         GI=GI+DI
0925         K3=K3+3
0926         IF ((ABS(CR)+ABS(DR)+ABS(CI)+ABS(DI)).GE.EPS) GOTO 71
0927         P=X*GR-Y*GI
0928         Q=X*GI+Y*GR
0929         GR=P
0930         GI=Q
0931         RETURN
0932         END
0933         SUBROUTINE SERAI(X,Y,AIR,AII)
0934 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0935 CCC   AIRY AI(Z), TAYLOR, COMPLEX Z      CCC
0936 CCC                                      CCC
0937 CCC   INPUTS:                            CCC
0938 CCC        X,Y,    REAL AND IMAGINARY    CCC
0939 CCC                PARTS OF Z            CCC
0940 CCC   OUTPUTS:                           CCC
0941 CCC        AIR,AII, REAL AND IMAGINARY   CCC
0942 CCC                 PARTS OF AI(Z)       CCC
0943 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0944         DOUBLE PRECISION X,Y,EPS,AIR,AII
0945         DOUBLE PRECISION FZR,FZI,GZR,GZI,CONS1,CONS2
0946         DOUBLE PRECISION D1MACH
0947         EPS=D1MACH(3)
0948         IF (EPS.LT.1.D-15) EPS=1.D-15
0949         CONS1=0.355028053887817239260D0
0950         CONS2=0.258819403792806798405D0
0951         CALL FG(X,Y,EPS,FZR,FZI,GZR,GZI)
0952         AIR=CONS1*FZR-CONS2*GZR
0953         AII=CONS1*FZI-CONS2*GZI
0954         RETURN
0955         END
0956         SUBROUTINE SERAID(X,Y,AIR,AII)
0957 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0958 CCC   AIRY AI'(Z), TAYLOR, COMPLEX Z     CCC
0959 CCC                                      CCC
0960 CCC   INPUTS:                            CCC
0961 CCC        X,Y,   REAL AND IMAGINARY     CCC
0962 CCC               PARTS OF Z             CCC
0963 CCC   OUTPUTS:                           CCC
0964 CCC        AIR,AII, REAL AND IMAGINARY   CCC
0965 CCC                 PARTS OF AI'(Z)      CCC
0966 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0967         DOUBLE PRECISION X,Y,EPS,AIR,AII
0968         DOUBLE PRECISION FZR,FZI,GZR,GZI,CONS1,CONS2
0969         DOUBLE PRECISION D1MACH
0970         EPS=D1MACH(3)
0971         IF (EPS.LT.1.D-15) EPS=1.D-15
0972         CONS1=0.355028053887817239260D0
0973         CONS2=0.258819403792806798405D0
0974         CALL FGP(X,Y,EPS,FZR,FZI,GZR,GZI)
0975         AIR=CONS1*FZR-CONS2*GZR
0976         AII=CONS1*FZI-CONS2*GZI
0977         RETURN
0978         END
0979         SUBROUTINE EXPAI(GAIR,GAII)
0980 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0981 CCC   AIRY AI(Z), ASYMPTOTIC EXPANSION, COMPLEX Z CCC
0982 CCC                                               CCC
0983 CCC   OUTPUTS:                                    CCC
0984 CCC        GAIR, GAII,  REAL AND IMAGINARY        CCC
0985 CCC                     PARTS OF AI(Z)            CCC
0986 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
0987         DOUBLE PRECISION EPS,GAIR,GAII
0988         DOUBLE PRECISION THET,R,TH15,S1,
0989      *  C1,R32,FACTO,TH025,S3,C3
0990         DOUBLE PRECISION DF1,PSIIR,PSIII,CK,DFRR,DFII,SUMAR,SUMAI,
0991      *  DFR,DFI,DELTAR,DELTAI,FAC1,FAC2
0992         DOUBLE PRECISION CO,DF
0993         DOUBLE PRECISION D1MACH
0994         INTEGER K
0995         DIMENSION CO(20)
0996         COMMON/PARAM3/THET,R,TH15,S1,C1,R32
0997         COMMON/PARAM4/FACTO,TH025,S3,C3
0998         SAVE CO
0999         DATA CO/-6.944444444444445D-2,3.713348765432099D-2,
1000      *  -3.799305912780064D-2,5.764919041266972D-2,-0.116099064025515D0,
1001      *  0.291591399230751D0,-0.877666969510017D0,3.07945303017317D0,
1002      *  -12.3415733323452D0,55.6227853659171D0,-278.465080777603D0,
1003      *  1533.16943201280D0,-9207.20659972641D0,59892.5135658791D0,
1004      *  -419524.875116551D0,3148257.41786683D0,-25198919.8716024D0,
1005      *  214288036.963680D0,-1929375549.18249D0,18335766937.8906D0/
1006         EPS=D1MACH(3)
1007         IF (EPS.LT.1.D-15) EPS=1.D-15
1008         DF1=1.5D0/R32
1009         PSIIR=C1
1010         PSIII=-S1
1011         K=0
1012         CK=1.D0
1013         DF=1.D0
1014         DFRR=1.D0
1015         DFII=0.D0
1016         SUMAR=1.D0
1017         SUMAI=0.D0
1018 80      DF=DF*DF1
1019         CK=CO(K+1)*DF
1020         DFR=DFRR
1021         DFI=DFII
1022         DFRR=DFR*PSIIR-DFI*PSIII
1023         DFII=DFR*PSIII+DFI*PSIIR
1024         DELTAR=DFRR*CK
1025         DELTAI=DFII*CK
1026         SUMAR=SUMAR+DELTAR
1027         SUMAI=SUMAI+DELTAI
1028         K=K+1
1029         IF (SUMAR.NE.0) THEN
1030           IF (ABS(DELTAR/SUMAR).GT.EPS)  GOTO 80
1031         ENDIF
1032         IF (SUMAI.NE.0) THEN
1033           IF (ABS(DELTAI/SUMAI).GT.EPS) GOTO 80
1034         ENDIF
1035         FAC1=FACTO*C3
1036         FAC2=FACTO*S3
1037         GAIR=FAC1*SUMAR+FAC2*SUMAI
1038         GAII=FAC1*SUMAI-FAC2*SUMAR
1039         RETURN
1040         END
1041         SUBROUTINE EXPAID(GAIR,GAII)
1042 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1043 CCC   AIRY AI'(Z), ASYMPTOTIC EXPANSION, COMPLEX Z CCC
1044 CCC                                                CCC
1045 CCC   OUTPUTS:                                     CCC
1046 CCC        GAIR, GAII,  REAL AND IMAGINARY         CCC
1047 CCC                     PARTS OF AI'(Z)            CCC
1048 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1049         DOUBLE PRECISION EPS,GAIR,GAII
1050         DOUBLE PRECISION THET,R,TH15,S1,
1051      *  C1,R32,FACTO,TH025,S3,C3
1052         DOUBLE PRECISION DF1,PSIIR,PSIII,VK,DFRR,DFII,SUMAR,SUMAI,
1053      *  DFR,DFI,DELTAR,DELTAI,FAC1,FAC2
1054         DOUBLE PRECISION CO,DF
1055         DOUBLE PRECISION D1MACH
1056         INTEGER K
1057         DIMENSION CO(20)
1058         COMMON/PARAM3/THET,R,TH15,S1,C1,R32
1059         COMMON/PARAM4/FACTO,TH025,S3,C3
1060         SAVE CO
1061         DATA CO/9.722222222222222D-2,-4.388503086419753D-2,
1062      *  4.246283078989484D-2,-6.266216349203230D-2,
1063      *  0.124105896027275D0,-0.308253764901079D0,
1064      *  0.920479992412945D0,-3.21049358464862D0,
1065      *  12.8072930807356D0,-57.5083035139143D0,
1066      *  287.033237109221D0,-1576.35730333710D0,
1067      *  9446.35482309593D0,-61335.7066638521D0,
1068      *  428952.400400069D0,-3214536.52140086D0,
1069      *  25697908.3839113D0,-218293420.832160D0,
1070      *  1963523788.99103D0,-18643931088.1072D0/
1071         EPS=D1MACH(3)
1072         IF (EPS.LT.1.D-15) EPS=1.D-15
1073         DF1=1.5D0/R32
1074         PSIIR=C1
1075         PSIII=-S1
1076         K=0
1077         DF=1.D0
1078         DFRR=1.D0
1079         DFII=0.D0
1080         SUMAR=1.D0
1081         SUMAI=0.D0
1082  81     DF=DF*DF1
1083         VK=CO(K+1)*DF
1084         DFR=DFRR
1085         DFI=DFII
1086         DFRR=DFR*PSIIR-DFI*PSIII
1087         DFII=DFR*PSIII+DFI*PSIIR
1088         DELTAR=DFRR*VK
1089         DELTAI=DFII*VK
1090         SUMAR=SUMAR+DELTAR
1091         SUMAI=SUMAI+DELTAI
1092         K=K+1
1093         IF (SUMAR.NE.0) THEN
1094           IF (ABS(DELTAR/SUMAR).GT.EPS)  GOTO 81
1095         ENDIF
1096         IF (SUMAI.NE.0) THEN
1097           IF (ABS(DELTAI/SUMAI).GT.EPS) GOTO 81
1098         ENDIF
1099         FAC1=FACTO*C3
1100         FAC2=FACTO*S3
1101         GAIR=-(FAC1*SUMAR-FAC2*SUMAI)
1102         GAII=-(FAC1*SUMAI+FAC2*SUMAR)
1103         RETURN
1104         END
1105
1106        SUBROUTINE BIZ(IFUN,IFAC,X0,Y0,GBIR,GBII,IERRO)
1107 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1108 C COMPUTATION OF THE AIRY FUNCTION BI(Z) OR ITS DERIVATIVE BI'(Z)
1109 C THE CODE USES THE CONNECTION OF BI(Z) WITH AI(Z).
1110 C                BIZ CALLS THE ROUTINE AIZ
1111 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1112 C  INPUTS:
1113 C    IFUN:
1114 C         * IFUN=1, THE CODE COMPUTES BI(Z)
1115 C         * IFUN=2, THE CODE COMPUTES BI'(Z)
1116 C    IFAC:
1117 C         * IFAC=1, THE CODE COMPUTES  BI(Z) OR BI'(Z)
1118 C         * IFAC=2, THE CODE COMPUTES NORMALIZED BI(Z) OR BI'(Z)
1119 C    X0:   REAL PART OF THE ARGUMENT Z
1120 C    Y0:   IMAGINARY PART OF THE ARGUMENT  Z
1121 C
1122 C  OUTPUTS:
1123 C    GBIR: REAL PART OF BI(Z) OR BI'(Z)
1124 C    GBII: IMAGINARY PART OF BI(Z) OR BI'(Z)
1125 C
1126 C    IERRO: ERROR FLAG
1127 C          * IERRO=0, SUCCESSFUL COMPUTATION
1128 C          * IERRO=1, COMPUTATION OUT OF RANGE
1129 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1130 C          MACHINE DEPENDENT CONSTANTS: FUNCTION D1MACH
1131 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1132 C   ACCURACY:
1133 C
1134 C     1) SCALED AIRY FUNCTIONS:
1135 C        RELATIVE ACCURACY BETTER THAN 10**(-13) EXCEPT CLOSE TO
1136 C        THE ZEROS, WHERE 10**(-13) IS THE ABSOLUTE PRECISION.
1137 C        GRADUAL LOSS OF PRECISION TAKES PLACE FOR |Z|>1000
1138 C        IN THE CASE OF PHASE(Z) CLOSE TO +3*PI/2 OR -3*PI/2.
1139 C     2) UNSCALED AIRY FUNCTIONS:
1140 C        THE FUNCTION OVERFLOWS/UNDERFLOWS FOR
1141 C        3/2*|Z|**(3/2)>LOG(OVER).
1142 C        FOR |Z|<30:
1143 C        A) RELATIVE ACCURACY FOR THE MODULUS (EXCEPT AT THE
1144 C           ZEROS) BETTER THAN 10**(-13).
1145 C        B) ABSOLUTE ACCURACY FOR MIN(R(Z),1/R(Z)) BETTER
1146 C           THAN 10**(-13), WHERE R(Z)=REAL(BI)/IMAG(BI)
1147 C           OR R(Z)=REAL(BI')/IMAG(BI').
1148 C        FOR |Z|>30, GRADUAL LOSS OF PRECISION TAKES PLACE
1149 C        AS |Z| INCREASES.
1150 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1151 C     AUTHORS:
1152 C        AMPARO GIL    (U. AUTONOMA DE MADRID, MADRID, SPAIN).
1153 C                      E-MAIL: AMPARO.GIL@UAM.ES
1154 C        JAVIER SEGURA (U. CARLOS III DE MADRID, MADRID, SPAIN).
1155 C                      E-MAIL: JSEGURA@MATH.UC3M.ES
1156 C        NICO M. TEMME (CWI, AMSTERDAM, THE NETHERLANDS).
1157 C                      E-MAIL: NICO.TEMME@CWI.NL
1158 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1159 C    REFERENCES:
1160 C         COMPUTING AIRY FUNCTIONS BY NUMERICAL QUADRATURE.
1161 C         NUMERICAL ALGORITHMS (2001).
1162 C         A. GIL, J. SEGURA, N.M. TEMME
1163 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1164        DOUBLE PRECISION X0,Y0,GBIR,GBII
1165        DOUBLE PRECISION OVER,UNDER,DL1,DL2,COVER,D1MACH
1166        DOUBLE PRECISION PI,PI3,PI23,SQRT3,C,S,C1,S1,U,V,X,Y,AR,AI,
1167      * APR,API,BR,BI,BPR,BPI,BBR,BBI,BBPR,BBPI,PHASE
1168        DOUBLE PRECISION THET,R,R32,THET32,A1,B1,DF1,EXPO,EXPOI,
1169      * ETAR,ETAI,ETAGR,ETAGI,PIHAL
1170        INTEGER IFUN,IFAC,IEXPF,IERR,IERRO
1171        COMMON/PARAM1/PI,PIHAL
1172        SQRT3=1.7320508075688772935D0
1173        PI=3.1415926535897932385D0
1174        PIHAL=1.5707963267948966192D0
1175        PI3=PI/3.D0
1176        PI23=2.D0*PI3
1177        X=X0
1178        C=0.5D0*SQRT3
1179        S=0.5D0
1180        IERRO=0
1181        IEXPF=0
1182        IF (Y0.LT.0.D0) THEN
1183          Y=-Y0
1184        ELSE
1185          Y=Y0
1186        ENDIF
1187        R=SQRT(X*X+Y*Y)
1188        R32=R*SQRT(R)
1189        THET=PHASE(X,Y)
1190        COVER=2.D0/3.D0*R32*ABS(COS(1.5D0*THET))
1191 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
1192        OVER=D1MACH(2)*1.D-3
1193 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
1194        UNDER=D1MACH(1)*1.D+3
1195        DL1=LOG(OVER)
1196        DL2=-LOG(UNDER)
1197        IF (DL1.GT.DL2) OVER=1/UNDER
1198        IF (IFAC.EQ.1) THEN
1199          IF (COVER.GE.LOG(OVER)) THEN
1200 CCC OVERFLOW/UNDERFLOW PROBLEMS.
1201 CCC   CALCULATION ABORTED
1202            IERRO=1
1203            GBIR=0
1204            GBII=0
1205          ENDIF
1206        ELSE
1207          IF (COVER.GE.(LOG(OVER)*0.2)) IEXPF=1
1208        ENDIF
1209        IF (IERRO.EQ.0) THEN
1210          IF (IFAC.EQ.1) THEN
1211            IF (Y.EQ.0.D0) THEN
1212 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1213 CCC  TAKE TWICE THE REAL PART OF EXP(-PI I/6) AI_(1)(Z)
1214 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1215              C1=-0.5D0
1216              S1=-0.5D0*SQRT3
1217              U=X*C1-Y*S1
1218              V=X*S1+Y*C1
1219              CALL AIZ(IFUN,IFAC,U,V,AR,AI,IERR)
1220              IF (IFUN.EQ.1) THEN
1221                BR=SQRT3*AR+AI
1222                BI=0.D0
1223              ELSE
1224                U=AR*C1-AI*S1
1225                V=AR*S1+AI*C1
1226                APR=U
1227                API=V
1228                BPR=SQRT3*APR+API
1229                BPI=0.D0
1230              ENDIF
1231            ELSE
1232              IF ((X.LT.0.D0).AND.(Y.LT.-X*SQRT3)) THEN
1233 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1234 CCC      2 PI/3  < PHASE(Z) < PI
1235 CCC      BI(Z)=EXP(I PI/6) AI_(-1)(Z) + EXP(-I PI/6) AI_(1)(Z)
1236 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1237                C1=-0.5D0
1238                S1=0.5D0*SQRT3
1239                U=X*C1-Y*S1
1240                V=X*S1+Y*C1
1241                CALL AIZ(IFUN,IFAC,U,V,AR,AI,IERR)
1242                IF (IFUN.EQ.1) THEN
1243                  BR=C*AR-S*AI
1244                  BI=C*AI+S*AR
1245                ELSE
1246                  U=AR*C1-AI*S1
1247                  V=AR*S1+AI*C1
1248                  APR=U
1249                  API=V
1250                  BPR=C*APR-S*API
1251                  BPI=C*API+S*APR
1252                ENDIF
1253 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1254 CCC   WE NEED ALSO AI_(1)(Z)
1255 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1256                C1=-0.5D0
1257                S1=-0.5D0*SQRT3
1258                U=X*C1-Y*S1
1259                V=X*S1+Y*C1
1260                CALL AIZ(IFUN,IFAC,U,V,AR,AI,IERR)
1261                IF (IFUN.EQ.1) THEN
1262                  S=-S
1263                  BR=BR+C*AR-S*AI
1264                  BI=BI+C*AI+S*AR
1265                ELSE
1266                  U=AR*C1-AI*S1
1267                  V=AR*S1+AI*C1
1268                  APR=U
1269                  API=V
1270                  S=-S
1271                  BPR=BPR+C*APR-S*API
1272                  BPI=BPI+C*API+S*APR
1273                ENDIF
1274              ELSE
1275 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1276 CCC   BI(Z) = I AI(Z) + 2 EXP(-I PI/6) AI_(1)(Z)
1277 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1278                C1=-0.5D0
1279                S1=-0.5D0*SQRT3
1280                U=X*C1-Y*S1
1281                V=X*S1+Y*C1
1282                CALL AIZ(IFUN,IFAC,U,V,AR,AI,IERR)
1283                IF (IFUN.EQ.1) THEN
1284                  BR=SQRT3*AR+AI
1285                  BI=-AR+SQRT3*AI
1286                ELSE
1287                  U=AR*C1-AI*S1
1288                  V=AR*S1+AI*C1
1289                  APR=U
1290                  API=V
1291                  BPR=SQRT3*APR+API
1292                  BPI=-APR+SQRT3*API
1293                ENDIF
1294                CALL AIZ(IFUN,IFAC,X,Y,AR,AI,IERR)
1295                IF (IFUN.EQ.1) THEN
1296                  BR=BR-AI
1297                  BI=BI+AR
1298                ELSE
1299                  BPR=BPR-AI
1300                  BPI=BPI+AR
1301                ENDIF
1302              ENDIF
1303            ENDIF
1304          ELSE
1305 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1306 CCC         SCALED BI AIRY FUNCTIONS         C
1307 CCC   WE USE THE FOLLOWING NORMALIZATION:    C
1308 CCC   LET ARGZ=ARG(Z), THEN:                 C
1309 CCC   A) IF  0 <= ARGZ <= PI/3               C
1310 CCC      BI=EXP(-2/3Z^3/2)BI                 C
1311 CCC   B) IF  PI/3 <= ARGZ <= PI              C
1312 CCC      BI=EXP(2/3Z^3/2)BI                  C
1313 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1314            THET=PHASE(X,Y)
1315            IF (THET.LE.PI3) THEN
1316              C1=-0.5D0
1317              S1=-0.5D0*SQRT3
1318              U=X*C1-Y*S1
1319              V=X*S1+Y*C1
1320              CALL AIZ(IFUN,IFAC,U,V,AR,AI,IERR)
1321              IF (IFUN.EQ.1) THEN
1322                BR=SQRT3*AR+AI
1323                BI=-AR+SQRT3*AI
1324              ELSE
1325                U=AR*C1-AI*S1
1326                V=AR*S1+AI*C1
1327                APR=U
1328                API=V
1329                BPR=SQRT3*APR+API
1330                BPI=-APR+SQRT3*API
1331              ENDIF
1332              IF (IEXPF.EQ.0) THEN
1333                R=SQRT(X*X+Y*Y)
1334                R32=R*SQRT(R)
1335                THET32=THET*1.5D0
1336                A1=COS(THET32)
1337                B1=SIN(THET32)
1338                DF1=4.D0/3.D0*R32
1339                EXPO=EXP(DF1*A1)
1340                EXPOI=1.D0/EXPO
1341                ETAR=EXPO*COS(DF1*B1)
1342                ETAI=EXPO*SIN(DF1*B1)
1343                ETAGR=EXPOI*COS(-DF1*B1)
1344                ETAGI=EXPOI*SIN(-DF1*B1)
1345                CALL AIZ(IFUN,IFAC,X,Y,AR,AI,IERR)
1346                IF (IFUN.EQ.1) THEN
1347                  BR=BR-AR*ETAGI-ETAGR*AI
1348                  BI=BI+AR*ETAGR-ETAGI*AI
1349                ELSE
1350                  BPR=BPR-AR*ETAGI-ETAGR*AI
1351                  BPI=BPI+AR*ETAGR-ETAGI*AI
1352                ENDIF
1353              ENDIF
1354            ENDIF
1355            IF ((THET.GT.PI3).AND.(THET.LE.PI23)) THEN
1356              IF (IEXPF.EQ.0) THEN
1357                R=SQRT(X*X+Y*Y)
1358                R32=R*SQRT(R)
1359                THET32=THET*1.5D0
1360                A1=COS(THET32)
1361                B1=SIN(THET32)
1362                DF1=4.D0/3.D0*R32
1363                EXPO=EXP(DF1*A1)
1364                EXPOI=1.D0/EXPO
1365                ETAR=EXPO*COS(DF1*B1)
1366                ETAI=EXPO*SIN(DF1*B1)
1367                ETAGR=EXPOI*COS(-DF1*B1)
1368                ETAGI=EXPOI*SIN(-DF1*B1)
1369                C1=-0.5D0
1370                S1=-0.5D0*SQRT3
1371                U=X*C1-Y*S1
1372                V=X*S1+Y*C1
1373                CALL AIZ(IFUN,IFAC,U,V,AR,AI,IERR)
1374                IF (IFUN.EQ.1) THEN
1375                  BBR=SQRT3*AR+AI
1376                  BBI=-AR+SQRT3*AI
1377                  BR=BBR*ETAR-BBI*ETAI
1378                  BI=BBI*ETAR+BBR*ETAI
1379                ELSE
1380                  U=AR*C1-AI*S1
1381                  V=AR*S1+AI*C1
1382                  APR=U
1383                  API=V
1384                  BBPR=SQRT3*APR+API
1385                  BBPI=-APR+SQRT3*API
1386                  BPR=BBPR*ETAR-BBPI*ETAI
1387                  BPI=BBPI*ETAR+BBPR*ETAI
1388                ENDIF
1389              ELSE
1390                IF (IFUN.EQ.1) THEN
1391                  BR=0.D0
1392                  BI=0.D0
1393                ELSE
1394                  BPR=0.D0
1395                  BPI=0.D0
1396                ENDIF
1397              ENDIF
1398              CALL AIZ(IFUN,IFAC,X,Y,AR,AI,IERR)
1399              IF (IFUN.EQ.1) THEN
1400                BR=BR-AI
1401                BI=BI+AR
1402              ELSE
1403                BPR=BPR-AI
1404                BPI=BPI+AR
1405              ENDIF
1406            ENDIF
1407            IF (THET.GT.PI23) THEN
1408              C1=-0.5D0
1409              S1=0.5D0*SQRT3
1410              U=X*C1-Y*S1
1411              V=X*S1+Y*C1
1412              CALL AIZ(IFUN,IFAC,U,V,AR,AI,IERR)
1413              IF (IFUN.EQ.1) THEN
1414                BR=C*AR-S*AI
1415                BI=C*AI+S*AR
1416              ELSE
1417                U=AR*C1-AI*S1
1418                V=AR*S1+AI*C1
1419                APR=U
1420                API=V
1421                BPR=C*APR-S*API
1422                BPI=C*API+S*APR
1423              ENDIF
1424              IF (IEXPF.EQ.0) THEN
1425 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1426 CCC   WE NEED ALSO AI_(1)(Z)
1427 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1428                 R=SQRT(X*X+Y*Y)
1429                 R32=R*SQRT(R)
1430                 THET32=THET*1.5D0
1431                 A1=COS(THET32)
1432                 B1=SIN(THET32)
1433                 DF1=4.D0/3.D0*R32
1434                 EXPO=EXP(DF1*A1)
1435                 EXPOI=1.D0/EXPO
1436                 ETAR=EXPO*COS(DF1*B1)
1437                 ETAI=EXPO*SIN(DF1*B1)
1438                 ETAGR=EXPOI*COS(-DF1*B1)
1439                 ETAGI=EXPOI*SIN(-DF1*B1)
1440                 C1=-0.5D0
1441                 S1=-0.5D0*SQRT3
1442                 U=X*C1-Y*S1
1443                 V=X*S1+Y*C1
1444                 CALL AIZ(IFUN,IFAC,U,V,AR,AI,IERR)
1445                 IF (IFUN.EQ.1) THEN
1446                   S=-S
1447                   BBR=C*AR-S*AI
1448                   BBI=C*AI+S*AR
1449                   BR=BR+ETAR*BBR-ETAI*BBI
1450                   BI=BI+BBI*ETAR+ETAI*BBR
1451                 ELSE
1452                   U=AR*C1-AI*S1
1453                   V=AR*S1+AI*C1
1454                   APR=U
1455                   API=V
1456                   S=-S
1457                   BBPR=C*APR-S*API
1458                   BBPI=C*API+S*APR
1459                   BPR=BPR+ETAR*BBPR-ETAI*BBPI
1460                   BPI=BPI+BBPI*ETAR+ETAI*BBPR
1461                 ENDIF
1462               ENDIF
1463             ENDIF
1464           ENDIF
1465           IF (Y0.LT.0) THEN
1466             BI=-BI
1467             BPI=-BPI
1468           ENDIF
1469           IF (IFUN.EQ.1) THEN
1470             GBIR=BR
1471             GBII=BI
1472           ELSE
1473             GBIR=BPR
1474             GBII=BPI
1475           ENDIF
1476         ENDIF
1477         RETURN
1478         END
1479
1480
1481
1482
1483
1484