[12] | 1 | subroutine alpham(dwtp,iend,np) |
---|
| 2 | c-------------------------------------------------------------------------- |
---|
| 3 | save |
---|
| 4 | c |
---|
| 5 | include 'param_sz.h' |
---|
| 6 | include 'constcom.h' |
---|
| 7 | include 'pcordcom.h' |
---|
| 8 | include 'syscom.h' |
---|
| 9 | include 'ucom.h' |
---|
| 10 | c |
---|
| 11 | common/aswap/acor(7,imaa) |
---|
| 12 | data itest/0/ |
---|
| 13 | c-------------------------------------------------------------------------- |
---|
| 14 | cbm !!! warning is for theta = 40.71 degres |
---|
| 15 | cbm data thetar,sr,cr,tr/.71052354,.65223071,.75802051,0.8604394/!LIU |
---|
| 16 | c-------------------------------------------------------------------------- |
---|
| 17 | c* |
---|
| 18 | c el(1,ne)=l distance travelled by the r.p. for a complete cycle |
---|
| 19 | c el(2,ne)=aper (vmax) for dropping |
---|
| 20 | c el(3,ne)=iout |
---|
| 21 | c el(4,ne)=wr energy of r.p. |
---|
| 22 | c el(5,ne)=thetd incident and exiting angles in degrees |
---|
| 23 | c el(6,ne)=wmax |
---|
| 24 | c el(7,ne)=deltab |
---|
| 25 | c el(8,ne)=l1 |
---|
| 26 | c el(9,ne)=l2 |
---|
| 27 | c el(10,ne)=xles |
---|
| 28 | c el(11,ne)=xhes |
---|
| 29 | c |
---|
| 30 | c ------------------------------- |
---|
| 31 | ne=rne |
---|
| 32 | if(itest.eq.0) then |
---|
| 33 | gr=el(4,ne)/erest |
---|
| 34 | bg=sqrt(gr*(2.+gr)) |
---|
| 35 | rhour=brhof*bg/el(7,ne) |
---|
| 36 | endif |
---|
| 37 | thetar=el(5,ne)*radian |
---|
| 38 | sr=sin(thetar) |
---|
| 39 | cr=cos(thetar) |
---|
| 40 | tr=tan(thetar) |
---|
| 41 | C --------------------------- |
---|
| 42 | if(gamma.le.1.5) go to 100 !added oct 20 89 |
---|
| 43 | c --------------------------- |
---|
| 44 | dphi=dcon*dwtp |
---|
| 45 | xp=bgx/bgz |
---|
| 46 | yp=bgy/bgz |
---|
| 47 | if(acor(7,np).eq.1.) go to 70 !particle in |
---|
| 48 | bz=bgz/gamma |
---|
| 49 | dz=bz*dphi |
---|
| 50 | x=x+xp*dz |
---|
| 51 | y=y+yp*dz |
---|
| 52 | z=z+dz |
---|
| 53 | if(acor(7,np).eq.-1.) then |
---|
| 54 | zaway=zloc(ne)-z !particle out |
---|
| 55 | zn=el(9,ne)-zaway |
---|
| 56 | u=-zn*cr-x*sr |
---|
| 57 | v=zn*sr-x*cr |
---|
| 58 | if(zaway.le.0.) then |
---|
| 59 | iend=1 |
---|
| 60 | acor(7,np)=0.0 |
---|
| 61 | endif |
---|
| 62 | return |
---|
| 63 | else |
---|
| 64 | zaway=zloc(ne-1)+el(8,ne)+x*tr-z |
---|
| 65 | zn=-zaway |
---|
| 66 | u=zn*cr-x*sr |
---|
| 67 | v=zn*sr+x*cr |
---|
| 68 | endif |
---|
| 69 | if(u.lt.0.) go to 65 |
---|
| 70 | bx=bgx/gamma |
---|
| 71 | by=bgy/gamma |
---|
| 72 | acor(1,np)=u |
---|
| 73 | acor(2,np)=v |
---|
| 74 | acor(3,np)=y |
---|
| 75 | acor(4,np)=bz*cr-bx*sr |
---|
| 76 | acor(5,np)=bz*sr+bx*cr |
---|
| 77 | acor(6,np)=by |
---|
| 78 | acor(7,np)=1. |
---|
| 79 | 65 continue |
---|
| 80 | return |
---|
| 81 | 70 continue |
---|
| 82 | c ----------------------- |
---|
| 83 | u=acor(1,np) |
---|
| 84 | v=acor(2,np) |
---|
| 85 | w=acor(3,np) |
---|
| 86 | bu=acor(4,np) |
---|
| 87 | bv=acor(5,np) |
---|
| 88 | bw=acor(6,np) |
---|
| 89 | c ------------------------------ |
---|
| 90 | if(abs(u).lt.0.1.and.(abs(v).ge.el(2,ne).or.abs(w).ge.el(6,ne))) |
---|
| 91 | . go to 100 |
---|
| 92 | c ------------------------------ |
---|
| 93 | theta1=atan(xp) |
---|
| 94 | c theta1=atan2(bx,bz) |
---|
| 95 | con=1.+xp**2+yp**2 |
---|
| 96 | c thetai=thetar+theta1 !injection angle for an electron |
---|
| 97 | c ------------------------- |
---|
| 98 | c buvsq=bu**2+bv**2 |
---|
| 99 | c buv=sqrt(buvsq) |
---|
| 100 | c rhou=bgd*buv !for an electron |
---|
| 101 | c |
---|
| 102 | c ac=deltab/(gamma*brhof) |
---|
| 103 | acd=dphi*el(7,ne)/(gamma*brhof) |
---|
| 104 | c phi0=0.5*(0.5*pi-thetai) |
---|
| 105 | c umax=2.*cos(phi0)*sqrt(buv/ac) |
---|
| 106 | c ------------------------------------ |
---|
| 107 | c h1=v*cr |
---|
| 108 | c hf=-h1-1.28*umax*theta1 |
---|
| 109 | c hf=-h1-1.28*umax*xp |
---|
| 110 | c |
---|
| 111 | dbu=acd*u*bv |
---|
| 112 | dbv=acd*(w*bw-u*bu) |
---|
| 113 | dbw=-acd*w*bv |
---|
| 114 | c |
---|
| 115 | bu=bu+dbu |
---|
| 116 | bv=bv+dbv |
---|
| 117 | bw=bw+dbw |
---|
| 118 | buvsq=bu**2+bv**2 |
---|
| 119 | buv=sqrt(buvsq) |
---|
| 120 | c |
---|
| 121 | du=bu*dphi |
---|
| 122 | dv=bv*dphi |
---|
| 123 | dw=bw*dphi |
---|
| 124 | u=u+du |
---|
| 125 | v=v+dv |
---|
| 126 | w=w+dw |
---|
| 127 | c print 110,x,z,u,v,bv |
---|
| 128 | c110 format(1x,5(e12.5,1x)) |
---|
| 129 | c 1.3=.85(r of d.t)/sr !!!!!!!! |
---|
| 130 | if(u.lt.1.3) go to 72 |
---|
| 131 | vsign=v/acor(2,np) |
---|
| 132 | if(vsign.lt.0.0.and.(u.le.el(10,ne) .or. u.ge.el(11,ne))) |
---|
| 133 | . go to 100 |
---|
| 134 | 72 continue |
---|
| 135 | rccc=(bu*dbv-bv*dbu)/dphi |
---|
| 136 | rhoa=abs(buvsq*buv/rccc) |
---|
| 137 | thetaa=atan2(bv,bu) |
---|
| 138 | c -------------------------- |
---|
| 139 | ur=u+x*sin(thetaa-theta1) !u of c.r. |
---|
| 140 | c if(ur.le.0.0) goto 66 |
---|
| 141 | rhor=rhour/ur |
---|
| 142 | rh=rhor+x |
---|
| 143 | ds=sqrt(du**2+dv**2) |
---|
| 144 | be=ds*cos(theta1) |
---|
| 145 | c be=ds !modified for test |
---|
| 146 | ang=be/rh |
---|
| 147 | dz=rhor*ang |
---|
| 148 | c |
---|
| 149 | ant=ang+theta1 |
---|
| 150 | sant=sin(ant) |
---|
| 151 | cant=cos(ant) |
---|
| 152 | c ct1=cos(theta1) |
---|
| 153 | arg=sant-(rh/rhoa)*sin(ang) |
---|
| 154 | if(abs(arg).ge.1.) go to 100 |
---|
| 155 | theta2=asin(arg) |
---|
| 156 | ct2=cos(theta2) |
---|
| 157 | h2=-rhor+rh*cos(ang)+rhoa*(ct2-cant) |
---|
| 158 | c --------- |
---|
| 159 | zn=z+dz |
---|
| 160 | z = zn |
---|
| 161 | x=h2 |
---|
| 162 | c y=y+yp*rhoa*(ant-theta2)*ct1 |
---|
| 163 | c yp=yp*ct1/ct2 |
---|
| 164 | y=w |
---|
| 165 | by=bw |
---|
| 166 | yp=by/(buv*ct2) |
---|
| 167 | xp=tan(theta2) |
---|
| 168 | bgz=bgz*sqrt(con/(1.+xp**2+yp**2)) |
---|
| 169 | bgx=xp*bgz |
---|
| 170 | bgy=yp*bgz |
---|
| 171 | gamma=sqrt(1.+bgx**2+bgy**2+bgz**2) |
---|
| 172 | c |
---|
| 173 | if(u.gt.0.) go to 90 |
---|
| 174 | acor(7,np)=-1. |
---|
| 175 | return |
---|
| 176 | 90 continue |
---|
| 177 | acor(1,np)=u |
---|
| 178 | acor(2,np)=v |
---|
| 179 | acor(3,np)=w |
---|
| 180 | acor(4,np)=bu |
---|
| 181 | acor(5,np)=bv |
---|
| 182 | acor(6,np)=bw |
---|
| 183 | return |
---|
| 184 | 100 continue |
---|
| 185 | c WRITE(NDIAG,101) IEND,(ACOR(KK,NP),KK=1,7) |
---|
| 186 | c write(ndiag,101) np,u,v,w,bu,bv,bw,acor(7,np) |
---|
| 187 | c write(ndiag,101) ne,x,y,z,bx,by,bz,gamma |
---|
| 188 | c101 format(1x,i5,7(e9.2,1x)) |
---|
| 189 | x=2.*el(2,ne) |
---|
| 190 | return |
---|
| 191 | end |
---|
| 192 | c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* |
---|