source: PSPA/parmelaPSPA/trunk/alpham.f @ 315

      subroutine alpham(dwtp,iend,np)   
c--------------------------------------------------------------------------
      save
c
      include 'param_sz.h'
      include 'constcom.h'
      include 'pcordcom.h'
      include 'syscom.h'
      include 'ucom.h'
c
      common/aswap/acor(7,imaa) 
      data itest/0/             
c--------------------------------------------------------------------------
cbm !!! warning is for theta = 40.71 degres
cbm      data thetar,sr,cr,tr/.71052354,.65223071,.75802051,0.8604394/!LIU 
c--------------------------------------------------------------------------
c*
c       el(1,ne)=l distance travelled by the r.p. for a complete cycle
c       el(2,ne)=aper (vmax) for dropping
c       el(3,ne)=iout 
c       el(4,ne)=wr energy of r.p.
c       el(5,ne)=thetd incident and exiting angles in degrees
c       el(6,ne)=wmax
c       el(7,ne)=deltab
c       el(8,ne)=l1
c       el(9,ne)=l2
c       el(10,ne)=xles
c       el(11,ne)=xhes
c
c       -------------------------------
      ne=rne
        if(itest.eq.0) then
                gr=el(4,ne)/erest
                bg=sqrt(gr*(2.+gr))
                rhour=brhof*bg/el(7,ne)
        endif           
        thetar=el(5,ne)*radian
        sr=sin(thetar)
        cr=cos(thetar)
        tr=tan(thetar)
C       ---------------------------
        if(gamma.le.1.5) go to 100      !added oct 20 89
c       ---------------------------
        dphi=dcon*dwtp
        xp=bgx/bgz
        yp=bgy/bgz
        if(acor(7,np).eq.1.) go to 70  !particle  in
        bz=bgz/gamma
        dz=bz*dphi
        x=x+xp*dz
        y=y+yp*dz
        z=z+dz
        if(acor(7,np).eq.-1.) then
                zaway=zloc(ne)-z        !particle out
                zn=el(9,ne)-zaway
                u=-zn*cr-x*sr
                v=zn*sr-x*cr
                if(zaway.le.0.) then
                        iend=1
                        acor(7,np)=0.0
                endif
                return
        else
                zaway=zloc(ne-1)+el(8,ne)+x*tr-z
                zn=-zaway
                u=zn*cr-x*sr
                v=zn*sr+x*cr
        endif
        if(u.lt.0.) go to 65
        bx=bgx/gamma
        by=bgy/gamma
        acor(1,np)=u
        acor(2,np)=v
        acor(3,np)=y
        acor(4,np)=bz*cr-bx*sr
        acor(5,np)=bz*sr+bx*cr
        acor(6,np)=by
        acor(7,np)=1.
65      continue
        return
70      continue
c       -----------------------
        u=acor(1,np)
        v=acor(2,np)
        w=acor(3,np)
        bu=acor(4,np)
        bv=acor(5,np)
        bw=acor(6,np)
c       ------------------------------
        if(abs(u).lt.0.1.and.(abs(v).ge.el(2,ne).or.abs(w).ge.el(6,ne)))
     .   go to 100
c       ------------------------------
        theta1=atan(xp)
c       theta1=atan2(bx,bz)
        con=1.+xp**2+yp**2
c       thetai=thetar+theta1    !injection angle for an electron
c       -------------------------
c       buvsq=bu**2+bv**2
c       buv=sqrt(buvsq)
c       rhou=bgd*buv    !for an electron
c
c       ac=deltab/(gamma*brhof)
        acd=dphi*el(7,ne)/(gamma*brhof)
c       phi0=0.5*(0.5*pi-thetai)
c       umax=2.*cos(phi0)*sqrt(buv/ac)
c       ------------------------------------
c       h1=v*cr
c       hf=-h1-1.28*umax*theta1
c       hf=-h1-1.28*umax*xp
c       
        dbu=acd*u*bv
        dbv=acd*(w*bw-u*bu)
        dbw=-acd*w*bv
c
        bu=bu+dbu
        bv=bv+dbv
        bw=bw+dbw
        buvsq=bu**2+bv**2
        buv=sqrt(buvsq)
c
        du=bu*dphi
        dv=bv*dphi
        dw=bw*dphi
        u=u+du
        v=v+dv
        w=w+dw
c       print 110,x,z,u,v,bv
c110    format(1x,5(e12.5,1x))
c       1.3=.85(r of d.t)/sr   !!!!!!!!
        if(u.lt.1.3) go to 72
        vsign=v/acor(2,np)
        if(vsign.lt.0.0.and.(u.le.el(10,ne) .or. u.ge.el(11,ne))) 
     .    go to 100
  72    continue
        rccc=(bu*dbv-bv*dbu)/dphi
        rhoa=abs(buvsq*buv/rccc)
        thetaa=atan2(bv,bu)
c       --------------------------
        ur=u+x*sin(thetaa-theta1)  !u of c.r.
c       if(ur.le.0.0) goto 66
        rhor=rhour/ur
        rh=rhor+x
        ds=sqrt(du**2+dv**2)
        be=ds*cos(theta1)
c       be=ds  !modified for test
        ang=be/rh
        dz=rhor*ang
c
        ant=ang+theta1
        sant=sin(ant)
        cant=cos(ant)
c       ct1=cos(theta1)
        arg=sant-(rh/rhoa)*sin(ang)
        if(abs(arg).ge.1.) go to 100
        theta2=asin(arg)
        ct2=cos(theta2)
        h2=-rhor+rh*cos(ang)+rhoa*(ct2-cant)
c       ---------
        zn=z+dz
        z = zn
        x=h2
c       y=y+yp*rhoa*(ant-theta2)*ct1
c       yp=yp*ct1/ct2
        y=w
        by=bw
        yp=by/(buv*ct2)
        xp=tan(theta2)
        bgz=bgz*sqrt(con/(1.+xp**2+yp**2))
        bgx=xp*bgz
        bgy=yp*bgz
        gamma=sqrt(1.+bgx**2+bgy**2+bgz**2)
c
        if(u.gt.0.) go to 90
        acor(7,np)=-1.
        return
  90    continue
        acor(1,np)=u
        acor(2,np)=v
        acor(3,np)=w
        acor(4,np)=bu
        acor(5,np)=bv
        acor(6,np)=bw
        return
 100    continue
c       WRITE(NDIAG,101) IEND,(ACOR(KK,NP),KK=1,7)
c       write(ndiag,101) np,u,v,w,bu,bv,bw,acor(7,np)
c       write(ndiag,101) ne,x,y,z,bx,by,bz,gamma
c101    format(1x,i5,7(e9.2,1x))
        x=2.*el(2,ne)
        return
        end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*