2 SUBROUTINE aiz(IFUN,IFAC,X0,Y0,GAIR,GAII,IERRO)
64 DOUBLE PRECISION x0,y0,gair,gaii,
x,
w,xd,wd
65 DOUBLE PRECISION over,under,dl1,dl2,cover,
d1mach
66 DOUBLE PRECISION pi,pihal,pih3,pisr,a,alf,thet,
r,th15,
s1,c1,
67 * r32,facto,th025,
s3,c3,f23,pi23,sqrt3,xa,ya,f23r,df1,df2,
68 * s11,c11,dex,dre,dima,gar,gai,
c,
s,u,
v,v0,ar,ai,ar1,ai1,
69 * ro,coe1,coe2,rex,dfr,dfi,ar11,ai11,
phase
70 INTEGER ifun,ifac,ierro,iexpf,iexpf2,
n
72 dimension xd(25),wd(25)
73 common/param1/
pi,pihal
74 common/param2/pih3,pisr,a,alf
75 common/param3/thet,
r,th15,
s1,c1,r32
76 common/param4/facto,th025,
s3,c3
79 DATA x,
w/.283891417994567679d-1,.170985378860034935d0,
80 *.435871678341770460d0,.823518257913030858d0,1.33452543254227372d0,
81 *1.96968293206435071d0,2.72998134002859938d0,3.61662161916100897d0,
82 *4.63102611052654146d0,5.77485171830547694d0,7.05000568630218682d0,
83 *8.45866437513237792d0,10.0032955242749393d0,11.6866845947722423d0,
84 *13.5119659344693551d0,15.4826596959377140d0,17.6027156808069112d0,
85 *19.8765656022785451d0,22.3091856773962780d0,24.9061720212974207d0,
86 *27.6738320739497190d0,30.6192963295084111d0,33.7506560850239946d0,
87 *37.0771349708391198d0,40.6093049694341322d0,.143720408803313866d0,
88 *.230407559241880881d0,.242253045521327626d0,.203636639103440807d0,
89 *.143760630622921410d0,.869128834706078120d-1,.4541750018329
90 * 15883d-1,.206118031206069497d-1,.814278821268606972d-2,.280266
91 *075663377634d-2,.840337441621719716d-3,.219303732907765020d-3,
92 *.497401659009257760d-4,.978508095920717661d-5,.166542824603725
93 *563d-5,.244502736801316287d-6,.308537034236207072d-7,.3332960
94 *72940112245d-8,.306781892316295828d-9,.239331309885375719d-10,
95 *.157294707710054952d-11,.864936011664392267d-13,.394819815
96 *638647111d-14,.148271173082850884d-15,.453390377327054458d-17/
97 DATA xd,wd/.435079659953445d-1,.205779160144678d0,
98 *.489916161318751d0,.896390483211727d0,1.42582496737580d0,
99 *2.07903190767599d0,2.85702335104978d0,3.76102058198275d0,
100 *4.79246521225895d0,5.95303247470003d0,7.24464710774066d0,
101 *8.66950223642504d0,10.2300817341775d0,11.9291866622602d0,
102 *13.7699665302828d0,15.7559563095946d0,17.8911203751898d0,
103 *20.1799048700978d0,22.6273004064466d0,25.2389175786164d0,
104 *28.0210785229929d0,30.9809287996116d0,34.1265753192057d0,
105 *37.4672580871163d0,41.0135664833476d0,.576354557898966d-1,
106 *.139560003272262d0,.187792315011311d0,.187446935256946d0,
107 *.150716717316301d0,.101069904453380d0,.575274105486025d-1,
108 *.280625783448681d-1,.117972164134041d-1,.428701743297432d-2,
109 *.134857915232883d-2,.367337337105948d-3,.865882267841931d-4,
110 *.176391622890609d-4,.309929190938078d-5,.468479653648208d-6,
111 *.607273267228907d-7,.672514812555074d-8,.633469931761606d-9,
112 *.504938861248542d-10,.338602527895834d-11,.189738532450555d-12,
113 *.881618802142698d-14,.336676636121976d-15,.104594827170761d-16/
115 pi=3.1415926535897932385d0
116 pihal=1.5707963267948966192d0
117 pih3=4.71238898038469d0
118 f23=.6666666666666666d0
119 pi23=2.09439510239320d0
120 pisr=1.77245385090552d0
121 sqrt3=1.7320508075688772935d0
128 IF (ya.LT.0.d0) ya=-ya
132 cover=2.d0/3.d0*r32*abs(cos(1.5d0*thet))
139 IF (dl1.GT.dl2) over=1/under
141 IF (cover.GE.
log(over))
THEN
148 IF (cover.GE.(
log(over)*0.2)) iexpf2=1
150 IF (cover.GE.(
log(over)*0.2)) iexpf=1
158 IF ((ya.LT.3.d0).AND.(xa.LT.1.3d0).AND.(xa.GT.-2.5d0))
THEN
160 CALL
serai(xa,ya,gar,gai)
174 gair=dre*gar-dima*gai
175 gaii=dre*gai+dima*gar
179 IF (y0.EQ.0.) gaii=0.d0
185 facto=0.5d0/pisr*
r**(-0.25d0)
186 IF (thet.GT.pi23)
THEN
204 IF (v0.LT.0.d0) ai1=-ai1
208 IF (iexpf2.EQ.0)
THEN
225 IF (v0.LT.0.d0) ai1=-ai1
230 ro=1.333333333333333d0*r32
238 gair=ar-(
c*ar11-
s*ai11)
239 gaii=ai-(
s*ar11+
c*ai11)
261 CALL
expai(gair,gaii)
265 a=0.1666666666666666d0
267 facto=0.280514117723058d0*
r**(-0.25d0)
269 IF (thet.LE.pihal)
THEN
278 IF ((thet.GT.pihal).AND.(thet.LE.pi23))
THEN
287 IF (thet.GT.pi23)
THEN
296 IF (thet.LE.pihal)
THEN
305 IF ((thet.GT.pihal).AND.(thet.LE.pi23))
THEN
314 IF (v0.LT.0.d0) ai1=-ai1
325 IF (thet.LE.pihal)
THEN
334 IF ((thet.GT.pihal).AND.(thet.LE.pi23))
THEN
343 IF (v0.LT.0.d0) ai1=-ai1
348 ro=1.333333333333333d0*r32
356 gair=ar-(
c*ar11-
s*ai11)
357 gaii=ai-(
s*ar11+
c*ai11)
372 gar=dre*gair-dima*gaii
373 gai=dre*gaii+dima*gair
376 IF (y0.EQ.0.) gaii=0.d0
381 alf=0.1666666666666666d0
382 facto=-0.270898621247918d0*
r**0.25d0
384 IF ((ya.LT.3.d0).AND.(xa.LT.1.3d0).AND.(xa.GT.-2.5d0))
THEN
386 CALL
seraid(xa,ya,gar,gai)
400 gair=dre*gar-dima*gai
401 gaii=dre*gai+dima*gar
405 IF (y0.EQ.0.) gaii=0.d0
411 facto=0.5d0/pisr*
r**0.25d0
412 IF (thet.GT.pi23)
THEN
430 IF (v0.LT.0.d0) ai1=-ai1
434 IF (iexpf2.EQ.0)
THEN
451 IF (v0.LT.0.d0) ai1=-ai1
456 ro=1.333333333333333d0*r32
464 gair=ar-(
c*ar11+
s*ai11)
465 gaii=ai-(-
s*ar11+
c*ai11)
490 IF (thet.LE.pihal)
THEN
497 CALL
airy1d(xd,wd,gair,gaii)
499 IF ((thet.GT.pihal).AND.(thet.LE.pi23))
THEN
506 CALL
airy2d(xd,wd,gair,gaii)
508 IF (thet.GT.pi23)
THEN
517 IF (thet.LE.pihal)
THEN
524 CALL
airy1d(xd,wd,ar1,ai1)
526 IF ((thet.GT.pihal).AND.(thet.LE.pi23))
THEN
533 CALL
airy2d(xd,wd,ar1,ai1)
535 IF (v0.LT.0.d0) ai1=-ai1
546 IF (thet.LE.pihal)
THEN
553 CALL
airy1d(xd,wd,ar1,ai1)
555 IF ((thet.GT.pihal).AND.(thet.LE.pi23))
THEN
562 CALL
airy2d(xd,wd,ar1,ai1)
564 IF (v0.LT.0.d0) ai1=-ai1
569 ro=1.333333333333333d0*r32
577 gair=ar-(
c*ar11+
s*ai11)
578 gaii=ai-(-
s*ar11+
c*ai11)
593 gar=dre*gair-dima*gaii
594 gai=dre*gaii+dima*gair
597 IF (y0.EQ.0) gaii=0.d0
602 IF (y0.LT.0.d0) gaii=-gaii
616 DOUBLE PRECISION x,
w,gair,gaii
617 DOUBLE PRECISION pih3,pisr,a,alf,thet,
r,th15,
s1,c1,
618 * r32,facto,th025,
s3,c3,sumar,sumai,df1,df1c1,
phi,phi6,
619 *
s2,
c2,dmodu,dmodu2,funr,funi,fac1,fac2,
phase
621 dimension
x(25),
w(25)
622 common/param2/pih3,pisr,a,alf
623 common/param3/thet,
r,th15,
s1,c1,r32
624 common/param4/facto,th025,
s3,c3
634 dmodu=sqrt(4.d0+df1*df1+4.d0*df1c1)
638 sumar=sumar+
w(
i)*funr
639 sumai=sumai+
w(
i)*funi
643 gair=fac1*sumar+fac2*sumai
644 gaii=fac1*sumai-fac2*sumar
658 DOUBLE PRECISION x,
w,gair,gaii
659 DOUBLE PRECISION pih3,pisr,a,alf,thet,
r,th15,
s1,c1,
660 * r32,facto,th025,
s3,c3,sumar,sumai,df1,df1c1,
phi,phi6,
661 *
s2,
c2,dmodu,dmodu2,funr,funi,fac1,fac2,
phase
662 DOUBLE PRECISION sqr2,sqr2i,
tau,tgtau,b,ang,ctau,cfac,ct,st,
663 * sumr,sumi,ttau,
beta
665 dimension
x(25),
w(25)
666 common/param2/pih3,pisr,a,alf
667 common/param3/thet,
r,th15,
s1,c1,r32
668 common/param4/facto,th025,
s3,c3
669 sqr2=1.41421356237310d0
670 sqr2i=0.707106781186548d0
682 df1=3.d0*
x(
i)/(ctau*r32)
683 df1c1=df1*sqr2i*0.5d0
690 dmodu=sqrt(4.d0+df1*df1*0.25d0-sqr2*df1)
697 sumar=cfac*(ct*sumr-st*sumi)
698 sumai=cfac*(ct*sumi+st*sumr)
701 gair=fac1*sumar+fac2*sumai
702 gaii=fac1*sumai-fac2*sumar
716 DOUBLE PRECISION x,
w,gair,gaii
717 DOUBLE PRECISION pih3,pisr,a,alf,thet,
r,th15,
s1,c1,
718 * r32,facto,th025,
s3,c3,sumar,sumai,df1,df1c1,
phi,phi6,
719 *
s2,
c2,dmodu,dmodu2,funr,funi,fac1,fac2,
phase
721 dimension
x(25),
w(25)
722 common/param2/pih3,pisr,a,alf
723 common/param3/thet,
r,th15,
s1,c1,r32
724 common/param4/facto,th025,
s3,c3
734 dmodu=sqrt(4.d0+df1*df1+4.d0*df1c1)
738 sumar=sumar+
w(
i)*funr
739 sumai=sumai+
w(
i)*funi
743 gair=fac1*sumar-fac2*sumai
744 gaii=fac1*sumai+fac2*sumar
758 DOUBLE PRECISION x,
w,gair,gaii
759 DOUBLE PRECISION pih3,pisr,a,alf,thet,
r,th15,
s1,c1,
760 * r32,facto,th025,
s3,c3,sumar,sumai,df1,df1c1,
phi,phi6,
761 *
s2,
c2,dmodu,dmodu2,funr,funi,fac1,fac2,
phase
762 DOUBLE PRECISION sqr2,sqr2i,
tau,tgtau,b,ang,ctau,cfac,ct,st,
763 * sumr,sumi,ttau,
beta
765 dimension
x(25),
w(25)
766 common/param2/pih3,pisr,a,alf
767 common/param3/thet,
r,th15,
s1,c1,r32
768 common/param4/facto,th025,
s3,c3
769 sqr2=1.41421356237310d0
770 sqr2i=0.707106781186548d0
782 df1=3.d0*
x(
i)/(ctau*r32)
783 df1c1=df1*sqr2i*0.5d0
790 dmodu=sqrt(4.d0+df1*df1*0.25d0-sqr2*df1)
797 sumar=cfac*(ct*sumr-st*sumi)
798 sumai=cfac*(ct*sumi+st*sumr)
801 gair=fac1*sumar-fac2*sumai
802 gaii=fac1*sumai+fac2*sumar
805 DOUBLE PRECISION FUNCTION phase(X,Y)
806 DOUBLE PRECISION pi,pihal,
x,
y,ay,
p
810 common/param1/
pi,pihal
811 IF ((
x.EQ.0).AND.(
y.EQ.0))
THEN
817 ELSEIF ((
x+ay).GE.0.d0)
THEN
826 SUBROUTINE fgp(X,Y,EPS,FR,FI,GR,GI)
840 DOUBLE PRECISION x,
y,
eps,fr,fi,
gr,gi
842 DOUBLE PRECISION x2,y2,u,
v,
p,
q,cr,
ci,dr,di
871 IF ((abs(cr)+abs(dr)+abs(
ci)+abs(di)).GE.
eps) goto 70
880 SUBROUTINE fg(X,Y,EPS,FR,FI,GR,GI)
896 DOUBLE PRECISION x2,y2,u,
v,
p,
q,cr,
ci,dr,di
897 DOUBLE PRECISION x,
y,
eps,fr,fi,
gr,gi
926 IF ((abs(cr)+abs(dr)+abs(
ci)+abs(di)).GE.
eps) goto 71
944 DOUBLE PRECISION x,
y,
eps,air,aii
945 DOUBLE PRECISION fzr,fzi,gzr,gzi,cons1,cons2
948 IF (
eps.LT.1.d-15)
eps=1.d-15
949 cons1=0.355028053887817239260d0
950 cons2=0.258819403792806798405d0
951 CALL
fg(
x,
y,
eps,fzr,fzi,gzr,gzi)
952 air=cons1*fzr-cons2*gzr
953 aii=cons1*fzi-cons2*gzi
967 DOUBLE PRECISION x,
y,
eps,air,aii
968 DOUBLE PRECISION fzr,fzi,gzr,gzi,cons1,cons2
971 IF (
eps.LT.1.d-15)
eps=1.d-15
972 cons1=0.355028053887817239260d0
973 cons2=0.258819403792806798405d0
975 air=cons1*fzr-cons2*gzr
976 aii=cons1*fzi-cons2*gzi
987 DOUBLE PRECISION eps,gair,gaii
988 DOUBLE PRECISION thet,
r,th15,
s1,
989 * c1,r32,facto,th025,
s3,c3
990 DOUBLE PRECISION df1,psiir,psiii,ck,dfrr,dfii,sumar,sumai,
991 * dfr,dfi,deltar,deltai,fac1,fac2
992 DOUBLE PRECISION co,df
996 common/param3/thet,
r,th15,
s1,c1,r32
997 common/param4/facto,th025,
s3,c3
999 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/
1007 IF (
eps.LT.1.d-15)
eps=1.d-15
1022 dfrr=dfr*psiir-dfi*psiii
1023 dfii=dfr*psiii+dfi*psiir
1029 IF (sumar.NE.0)
THEN
1030 IF (abs(deltar/sumar).GT.
eps) goto 80
1032 IF (sumai.NE.0)
THEN
1033 IF (abs(deltai/sumai).GT.
eps) goto 80
1037 gair=fac1*sumar+fac2*sumai
1038 gaii=fac1*sumai-fac2*sumar
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
1058 common/param3/thet,
r,th15,
s1,c1,r32
1059 common/param4/facto,th025,
s3,c3
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/
1072 IF (
eps.LT.1.d-15)
eps=1.d-15
1086 dfrr=dfr*psiir-dfi*psiii
1087 dfii=dfr*psiii+dfi*psiir
1093 IF (sumar.NE.0)
THEN
1094 IF (abs(deltar/sumar).GT.
eps) goto 81
1096 IF (sumai.NE.0)
THEN
1097 IF (abs(deltai/sumai).GT.
eps) goto 81
1101 gair=-(fac1*sumar-fac2*sumai)
1102 gaii=-(fac1*sumai+fac2*sumar)
1106 SUBROUTINE biz(IFUN,IFAC,X0,Y0,GBIR,GBII,IERRO)
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
1182 IF (y0.LT.0.d0)
THEN
1190 cover=2.d0/3.d0*r32*abs(cos(1.5d0*thet))
1197 IF (dl1.GT.dl2) over=1/under
1199 IF (cover.GE.
log(over))
THEN
1207 IF (cover.GE.(
log(over)*0.2)) iexpf=1
1209 IF (ierro.EQ.0)
THEN
1219 CALL
aiz(ifun,ifac,u,
v,ar,ai,ierr)
1232 IF ((
x.LT.0.d0).AND.(
y.LT.-
x*sqrt3))
THEN
1241 CALL
aiz(ifun,ifac,u,
v,ar,ai,ierr)
1260 CALL
aiz(ifun,ifac,u,
v,ar,ai,ierr)
1282 CALL
aiz(ifun,ifac,u,
v,ar,ai,ierr)
1294 CALL
aiz(ifun,ifac,
x,
y,ar,ai,ierr)
1315 IF (thet.LE.pi3)
THEN
1320 CALL
aiz(ifun,ifac,u,
v,ar,ai,ierr)
1332 IF (iexpf.EQ.0)
THEN
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)
1347 br=br-ar*etagi-etagr*ai
1348 bi=bi+ar*etagr-etagi*ai
1350 bpr=bpr-ar*etagi-etagr*ai
1351 bpi=bpi+ar*etagr-etagi*ai
1355 IF ((thet.GT.pi3).AND.(thet.LE.pi23))
THEN
1356 IF (iexpf.EQ.0)
THEN
1365 etar=expo*cos(df1*
b1)
1366 etai=expo*sin(df1*
b1)
1367 etagr=expoi*cos(-df1*
b1)
1368 etagi=expoi*sin(-df1*
b1)
1373 CALL
aiz(ifun,ifac,u,
v,ar,ai,ierr)
1377 br=bbr*etar-bbi*etai
1378 bi=bbi*etar+bbr*etai
1386 bpr=bbpr*etar-bbpi*etai
1387 bpi=bbpi*etar+bbpr*etai
1398 CALL
aiz(ifun,ifac,
x,
y,ar,ai,ierr)
1407 IF (thet.GT.pi23)
THEN
1412 CALL
aiz(ifun,ifac,u,
v,ar,ai,ierr)
1424 IF (iexpf.EQ.0)
THEN
1436 etar=expo*cos(df1*
b1)
1437 etai=expo*sin(df1*
b1)
1438 etagr=expoi*cos(-df1*
b1)
1439 etagi=expoi*sin(-df1*
b1)
1444 CALL
aiz(ifun,ifac,u,
v,ar,ai,ierr)
1449 br=br+etar*bbr-etai*bbi
1450 bi=bi+bbi*etar+etai*bbr
1459 bpr=bpr+etar*bbpr-etai*bbpi
1460 bpi=bpi+bbpi*etar+etai*bbpr