program wibz
	complex ax,f1,f0,y1,y2,ya,yb,z,dic,zis
        dimension am(21),fi(21),t1(18),fa(18),fb(18) 
	dimension ax(143),is(21)
        common/bl1/d0,f,b0,dic 
        common/bl2/r0,al,ds,ha 
        common/bl3/pi,h1,e1,h,hd,e2,h2,zis 
        common/bl4/n,nl,nt,na
        external fsl 
        hd=0.18
        e1=9.6
        e2=9.6
        h2=0.18
        h=0.44
        b0=1.57
        ds=0.15
        r0=0.1
        d0=0.088
        ha=1.0
        nt=8
        na=10
        n=16
        al=2.1
        zis=cmplx(0.0,0.2)
        pi=3.141592
        h1=2*al/n
        x=sqrt(abs(ds*ds+d0*d0))
	di=alog(4.)*alog((x+ds)/
     *	(x-ds))-2.*ds*alog(ds)/x+2.*ds*alog(x)/x
        dic=cmplx(4./(e1+1.)*di,ds/(e1+1.))
        call asg4c(ax,n+1,1,fsl,is)
        nz=n/2+1
        ya=((r0-h1)*ax(nz+2)-(r0-2.*h1)*ax(nz+1))/h1
        yb=cmplx(0.,1.)
        x=sf24r(r0)
        z=-(yb/ya*x)
        i1=n+1
        do 10 i=1,i1
        br=real(ax(i))
        bi=aimag(ax(i))
        am(i)=sqrt(abs(br*br+bi*bi))
10      fi(i)=atan(bi/br)
        ya=(0.,0.)
        yb=(0.,0.)
        dt=pi/(2.*nt)
        do 20 j=1,nt
        t1(j)=j*dt+0.01
        tj=t1(j)
        call ft(tj,f0,f1)
        do 30 i=2,n
        xi=-al+(i-1)*h1
        x1=cos(xi*sin(tj))
        x2=sin(xi*sin(tj))
        y1=cmplx(x1,x2)
        y2=cmplx(0.,1.)
        ya=ya+ax(i)*(y1*f0*cos(tj)+y2*(f1*sin(tj)))
30      yb=yb+ax(i)*f0
        ya=ya*cos(tj)*h1
        yb=yb*cos(tj)*h1
        fa(j)=cabs(ya)
        fb(j)=cabs(yb)
        open(1,file='fa.txt',access='append')
        write(1,*)fa(j)
        open(2,file='fb.txt',access='append')
        write(2,*)fb(j)
20      continue
        open(3,file='ax.txt')
        write(3,*)(ax(i),i=1,i1)
        open(4,file='am.txt')
        write(4,*)(am(i),i=1,i1)
        open(5,file='fi.txt')
        write(5,*)(fi(i),i=1,i1)
        open(16,file='z.txt',access='append')
        write(16,*)z
        open(17,file='t1.txt')
        write(17,*)t1
        open(18,file='h1.txt',access='append')
        write(18,*)h1
        stop
        end

	subroutine fsl(a,k)
        complex a,dic
        dimension a(1)
        common/bl1/d0,sj,b0,dic
        common/bl2/r0,al,ds,ha
        common/bl3/pi,h1,e1,h,hd,e2,h2,zis
        common/bl4/n,nl,nt,na
        external g1
	complex gaus3,y,y1,y2,zis
        sj=(k-1)*h1-al
        ah=1./2./h1/h1
        in1=n-1
        do 10 i=3,in1,2
        si=(i-1)*h1-al
        s1=si-2.*h1
        s2=si-h1
        s3=si+h1
        s4=si+2.*h1
        if (i-k) 11,12,11
11      y1=gaus3(s1,si,s1,s2,g1)
        y2=gaus3(si,s4,s3,s4,g1)
        y=ah*(y1+y2)
        a(i)=y/2.
        goto 10
12      y1=gaus3(s1,(si-ds),s1,s2,g1)
        y2=gaus3((si+ds),s4,s3,s4,g1)
        y=ah*(y1+y2)
        a(i)=y/2.+dic
10      continue
        do 20 i=2,n,2
        si=(i-1)*h1-al
        s2=si-h1
        s3=si+h1
        if (i-k) 21,22,21
21      y1=gaus3(s2,s3,s2,s3,g1)
        a(i)=y1*(-2.*ah)
        goto 20
22      y1=gaus3(s2,(si-ds),s2,s3,g1)
        y2=gaus3((si+ds),s3,s2,s3,g1)
        y=-2.*ah*(y1+y2)
        a(i)=y+dic
20      continue
        a(1)=-sin(sj)
        a(n+1)=-cos(sj)
        a(n+2)=-b0*sin(abs(sj))/60.
        return
        end

        subroutine fg(a1,fg1,fg2)
        common/bl3/pi,h1,e1,h,hd,e2,h2,zis
        complex y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,y11,y13,y14
        complex r0h,r1h,r0e,r1e,fg1,fg2,zis
        a=a1*a1
        if (a.ge.1) goto 1
        ci=sqrt(abs(1.-a))
        cr=0.
        goto 2
1       cr=sqrt(abs(a-1.))
        ci=0.
2       y0=cmplx(cr,ci)
        if (a.ge.e1) goto 11
        di=sqrt(abs(e1-a))
        dr=0.
        goto 12
11      dr=sqrt(abs(a-e1))
        di=0.
12      y1=cmplx(dr,di)
        if (a.ge.e2) goto 21
        bi=sqrt(abs(e2-a))
        br=0.
        goto 22
21      br=sqrt(abs(a-e2))
        bi=0.
22      y2=cmplx(br,bi)
        y13=cexp(-2.*y0*h)
        y14=cexp(-2.*y1*hd)
        y5=cexp(-2.*y2*h2)
        y3=(1.-y13)/(1.+y13)
        y4=(1.-y14)/(1.+y14)
        y6=y1*(zis-y3)-y0*(1.+zis*y3)*y4
        y7=y0*(1.+zis*y3)+y1*(zis+y3)*y4
        y8=e1*y0*(1.-zis*y3)-y1*(zis-y3)*y4
        y9=y1*(zis-y3)-e2*y0*(1.-zis*y3)*y4
        r0h=(y2*y6-y1*y7)/(y2*y6+y1*y7)
        r1h=(y2-y0)/(y2+y0)
        r0e=(e2*y1*y8+e1*y2*y9)/
     *  (e2*y1*y8-e1*y2*y9)
        r1e=(e2*y0-y2)/(e2*y0+y2)
        y10=(r0h*(1.+r1h*y5)+r1h*y5*(1.+r0h))/
     *  (1.-r0h*r1h*y5)*a/y0
	y11=-(1.+r0h+r1h*(1.+r0h)*y5)*(e2-e1)/
     *	(e2+e1)/(1.-r1h*r0h*y5)*a/y2-y2/
     *  a*(r0h-r1h*(1.+r0h)*y5)/
     *	(1.-r1h*r0h*y5)+y2/
     *  a*(r0e-r1e*(1.-r0e)*y5)/(1.-r1e*r0e*y5)
        fg1=y10
        fg2=y11
        return
        end

        complex function g1(s0)
        common/bl1/d0,sj,b0,dic
        common/bl2/r0,al,ds,ha
        common/bl3/pi,h1,e1,h,hd,e2,h2,zis
        common/bl4/n,nl,nt,na
	complex y1,y2,y3,y4,y5,dic,zis
	complex fg1,fg3,fg4,fg6,fr1
        a=0.57735
        y1=(0.,0.)
        y2=(0.,0.)
        y3=(0.,0.)
        y4=(0.,0.)
        y5=(0.,0.)
        r1=abs(sj-s0)
        if (r1.le.0.01) return
        call fr(sj,s0,fr1)
        do 1 j=1,na
        xj=(2*j-1)*ha/2.
        a1=xj-ha/2.*a
        a2=xj+ha/2.*a
        p1=sf24r(a1*r1)
        p3=sf24r(a2*r1)
        call fg(a1,fg1,fg3)
        call fg(a2,fg4,fg6)
        y1=y1+p1*fg1+p3*fg4
1       y5=y5+p1*fg3+p3*fg6
        g1=2./(e1+1.)*(fr1+y1*ha/2.)
     *  +(y3+y5)*ha/2.
        return
        end

        subroutine fr(s,s0,fr1)
        common/bl1/d0,f,b0,dic
        common/bl3/pi,h1,e1,h,hd,e2,h2,zis
	complex fr1,y1,y3,dic,zis
        r1=abs(s-s0)
        x1=sqrt(abs(r1*r1+d0*d0))
        x3=fe(d0/x1)
        y1=cmplx(cos(r1),-sin(r1))
        y3=cmplx(2./pi*x3,2./pi*r1*x3-x1)
        fr1=y1*y3/x1
        return
        end

        complex function gaus3(a,b,c,d,f)
        complex f
        b1=(b+a)/2.
        b2=(b-a)/2.
        x1=b1-b2*0.861
        x2=b1-b2*0.34
        x3=b1+b2*0.34
        x4=b1+b2*0.861
	gaus3=b2*((f(x1)*(x1-c)*(x1-d)+f(x4)*(x4-c)*(x4-d))*0.348
     *	+(f(x2)*(x2-c)*(x2-d)+f(x3)*(x3-c)*(x3-d))*0.652)
        return
        end

        function fe(a)
        a1=1.-a*a
        a2=a1*a1
        a3=a2*a2
        b1=1.38629
        c1=0.5
        b2=0.09657
        c2=0.125
        b3=0.03095
        c3=0.0703
        b4=0.01694
        c4=0.04816
        b5=0.01974
        c5=0.03072
        b6=0.01724
        c6=0.0105
        b7=0.00305
        c7=0.00079
	fe=b1+b2*a1+b3*a2+b4*a2*a1+b5*a3+b6*a3*a1+b7*a3*a2-alog(a1)*
     *	(c1+c2*a1+c3*a2+c4*a1*a2+c5*a3+c6*a3*a1+c7*a3*a2)
        return
        end

        subroutine besl0r(x,z)
        real x,t,z
        x=abs(x)
        if (x.le.4) goto 1
        t=16./x/x
        z=2.*(cos(x-0.7853981634)*(((((-0.0000037*t*
     *	0.0000173565)*t-0.0000487613)*t+0.00017343)*
     *	t-0.001753)*t+0.3989422793)-sin(x-0.78539816)*
     *	(((((0.0000032312*t-0.0000142078)*t+0.0000342468)*
     *	t-0.0000869791)*t+0.0004564324)*t-0.01246694)*4./x)
     *  /sqrt(x)
        goto 2
1       t=x*x/16
        z=((((((-0.00050144*t+0.0076771853)*t-0.07092535)*t*
     *	0.4443584263)*t-1.7777560599)*t+3.9999973)
     *  *t-3.9999998721)*t+1
2       continue
        return
        end

        function sf24r(x)
        call besl0r(x,z)
        sf24r=z
        return
        end

        subroutine asg4c(a,n,m,fsl,ir)
        complex a,am,r
        integer n,m,ir,l,i,j,k,nm,nx,ls,kn,ke,kl,ks
        integer n1,n2,n3,nj,nk,na,nz
        real r1,r2,cabs
        dimension a(1),ir(n) 
        nm=n+m
        if (m-n) 1,2,2
1       l=nm*nm/4+nm
        goto 3  
2       l=n*(m+2)-1
3       ke=0
        ls=0
4       ks=ke
        kn=ks+1
        kl=(l-ke*(nm-ke))/nm
        nx=ls+1
        ke=ks+kl
        if (n-ke) 5,6,6
5       ke=n
        kl=n-ks
6       nz=nm-ks
        do 9 k=kn,ke
        call fsl(a(nx),n,m,k)
        if (ks.eq.0) goto 9
        n1=nx
        n2=n-1
        do 7 i=1,ks
        j=ir(i)
        nj=n2+j
        r=a(n1)  
        a(n1)=a(nj)
        a(nj)=r
7       n1=n1+1
        nk=n2+nm
        na=0
        n3=n1-1
        do 8 i=nx,n3
        r=a(i) 
        do 8 j=n1,nk
        na=na+1
8       a(j)=a(j)-r*a(na)
9       nx=nx+nm
        n2=ls 
        do 24 k=kn,ke
        nx=n2+k
        n3=n2+nm
        r1=0.
        nk=k
        do 10 nj=k,n
        r2=cabs(a(nx))
        nx=nx+1
        if (r1.ge.r2) goto 10
        r1=r2
        nk=nj
10      continue
        ir(k)=nk
        nj=n2+nk
        am=1./a(nj)
        nx=n2+k
        do 11 na=nx,n3
11      a(na)=am*a(na)
        if (ks.eq.0) goto 25
        n1=k-ks
        na=nk-ks
        i=(ks-1)*nz+n1
        do 12 j=n1,i,nz
        r=a(j)
        a(j)=a(na)
        a(na)=r
12      na=na+nz
25      n1=ls+k
        na=ls+nk
        i=(kl-1)*nm+n1
        do 13 j=n1,i,nm
        r=a(j)
        a(j)=a(na)
        a(na)=r
13      na=na+nm
        n1=nx+1
        if (ks.eq.0) goto 26
        nj=k-ks
        i=(ks-1)*nz+nj
        do 14 j=nj,i,nz
        r=a(j)
        na=j
        do 14 nk=n1,n3
        na=na+1
14      a(na)=a(na)-r*a(nk)
26      nj=ls+k
        i=(kl-1)*nm+nj
        do 16 j=nj,i,nm
        if (j.eq.nx) goto 16
        r=a(j)
        na=j
        do 15 nk=n1,n3
        na=na+1
15      a(na)=a(na)-r*a(nk)
16      continue
24      n2=n2+nm
        n3=nm-ke
        n1=(ks-1)*n3+1
        if (ks.eq.0) goto 27
        nj=kl+1
        nk=n3
        do 18 k=1,n1,n3
        n2=nj
        do 17 i=k,nk
        a(i)=a(n2)
17      n2=n2+1
        nk=nk+n3
18      nj=nj+nz
27      n2=n1+kl*n3
        n1=n1+n3
        nj=ls+ke+1
        nk=n1+n3-1
        do 20 k=n1,n2,n3
        na=nj
        do 19 i=k,nk
        a(i)=a(na)
19      na=na+1
        nk=nk+n3
20      nj=nj+nm
        ls=ke*n3
        if (ke.ne.n) goto 4
        if (n.eq.1) goto 23
        n1=n-1
        n3=n1*m
        n2=n3-m+1
        do 22 k=1,n1
        kl=n-k
        j=ir(kl)
        nj=(j-1)*m
        do 21 i=n2,n3
        nj=nj+1
        r=a(i)
        a(i)=a(nj)
21      a(nj)=r
        n3=n2-1
22      n2=n2-m
23      return
        end

        subroutine ft(t,f0,f1)
        common/bl3/pi,h1,e1,h,hd,e2,h2,zis
        common/bl4/n,nl,nt,na
        complex y1,y2,y3,y4,y5,y6,y7,r0h,r0e,r1h,r1e,f0,f1,zis
        a=cos(t)
        b=sin(t)
        x0=a
        x1=sqrt(abs(e1-b*b))
        x2=sqrt(abs(e2-b*b))
        x3=x0*h
        x4=x1*hd
        x5=x2*h2
        x6=sin(x3)/cos(x3)
        x7=sin(x4)/cos(x4)
        y1=exp(-x5)
        x10=x1*x6+x0*x7
        x11=x0-x1*x6*x7
        y2=cmplx(x0*x10,-x1*x11)
        y3=cmplx(x2*x10,x1*x11)
        x12=x1*x6+e1*x0*x7
        x13=e1*x0-x1*x7*x6
        y4=cmplx(-e1*x2*x12,e2*x1*x12)
        y5=cmplx(e1*x2*x12,e2*x1*x12)
        r0h=y2/y3
        r0e=-y4/y5
        r1h=(x2-x0)/(x2+x0)
        r1e=(e2*x0-x2)/(e2*x0+x2)
        y6=1.-r1h*r0h*y1
        y7=1.-r1e*r0e*y1
        f0=(1.-r1h)/y6*(1.+r0h)
        f1=1./e2*(1.+r1e)/y7*(1.-r0e)*0.5
        return
        end