program angles, xst,scatflg,thp,php,j) !bec implicit none include 'const.inc' integer i,j,loop ! bec real vec1(1), vec(3), qrand(1), HRNDM1 double precision en0, x0, y0, z0, thx0, thy0, tx1, tx2 double precision ty1, ty2, tz, x, y, z, thx, thy, ep, enp double precision rfx, rfy, rfz double precision dxs, T, hedens, cs, qp, eqp, th_test double precision thsc, phsc, thxsc, thysc double precision henumb, dreact, domga double precision th0, ph0, domega, rdist(3) double precision k0, kp, thp, php, thxp, thyp, xdis, ydis double precision qp1, qp2, csthsc1, csthsc2 integer qmin, qmax double precision q, eq, wmin, wmax, wtest, xst, xstot double precision denat, rc, rg, dsig0, ratio, thxtst, thytst double precision c1,c2,c3,con1,con2,con3,mm, thp1, thp2 double precision norm, wp, sinphp, tharg, pharg logical scatflg c call hlimit(max_paw) c call rluxgo(3,120566405,0,0) ! RANLUX luxury level 2 call hbook1(1,'th_o (deg) ',200,0.,360.,0.) call hbook1(2,'phi_o (deg) ',200,0.,360.,0.) call hbook1(3,'th_s (deg) ',200,0.,360.,0.) call hbook1(4,'phi_s (deg) ',200,0.,360.,0.) do i = 1, 100 321 call ranlux(vect,4) thindex = 0.003*5.0d0 ! 5 Ang neutron if (abs(1.d0-(1.d0-dcos(thindex))*vect(1)).gt.1.0d0) goto 321 th = dacos(1.d0-(1.d0-dcos(thindex))*vect(1)) ! isotropic (0,thindex) bec ph = 2.0d0*pi*vect(2) ! azimuthal angle (radians) thx0 = datan(dtan(th)*dcos(ph)) thy0 = datan(dtan(th)*dsin(ph)) if(thx0.ge.0..and.thy0.ge.0.) then ph0 = datan(dtan(thy0)/dtan(thx0)) elseif(thx0.lt.0..and.thy0.ge.0.) then ph0 = pi + datan(dtan(thy0)/dtan(thx0)) elseif(thx0.lt.0..and.thy0.lt.0.) then ph0 = pi + datan(dtan(thy0)/dtan(thx0)) elseif(thx0.ge.0..and.thy0.lt.0.) then ph0 = 2.*pi + datan(dtan(thy0)/dtan(thx0)) endif th0 = datan(sqrt(dtan(thx0)**2 + dtan(thy0)**2)) write(6,*) ' th,ph,th0,ph0 ', th,ph,th0,ph0 cc here is stuff for k, kp, thsc if want to put it back in c call ranlux(vec,3) thsc = 2.0*pi*vec(1) phsc = 2.0*pi*vec(2) c cc php = datan(dtan(phsc)*dcos(thx0)/dcos(thy0)) cc above is wrong since by definition cos(thx0)=cos(thy0) cc also, the original stuff below produces phi values only from 0-90 cc php = phsc ! bec cc if (abs(dsin(thsc)*dcos(phsc)/(dcos(php)*dcos(thx0)*dcos(th0))) cc & .gt.1.0) then cc write(*,*) 'thsc:',thsc,' phsc:',phsc,' php:',php cc & ,' arg:',dsin(thsc)*dcos(phsc)/(dcos(php)*dcos(thx0)*dcos(th0)) cc thp = dasin(1.0d0) cc else cc thp = dasin(dsin(thsc)*dcos(phsc)/(dcos(php)*dcos(thx0)* cc & dcos(th0))) cc endif c c Using transformation matrix from x,y,z to X,Y,Z (lab) where k-vect is along c z-axis. Rotate by phi about Z, then theta about y, and chi=0 about z (no need for this rot). c This rotation creates CS, x,y,z relative to X,Y,Z, where theta and phi locate z axis c in X,Y,Z (lab) frame. We'll take k-vector along z, so c theta_o and phi_o equal theta and phi of rotated frame. kp (scattered vec) is c theta_sc from z (=k) with azimuthal angle phi_sc. Use following to rotate back to X,Y,Z, giving c theta_p and phi_p in X,Y,Z. c th0 = datan( dsqrt(dtan(thx0)**2+dtan(thy0)**2) ) ! from current thx, thy ph0 = datan(dtan(thy0)/dtan(thx0)) tharg = -dsin(th0)*dsin(thsc)*dcos(phsc)+dcos(th0)*dcos(thsc) if (tharg.gt.1.0d0) then write(*,*) ' theta_p error : cos(th_p) =',tharg tharg=1.0d0 else if (tharg.lt.-1.0d0) then write(*,*) ' theta_p error : cos(th_p) =',tharg tharg=-1.0d0 endif thp = dacos(tharg) pharg = (dsin(thsc)*dcos(phsc)*dcos(ph0)*dcos(th0) - & dsin(thsc)*dsin(phsc)*dsin(ph0) + & dcos(thsc)*dcos(ph0)*dsin(th0) )/dsin(thp) if (pharg.gt.1.0d0) then write(*,*) ' phi_p error : cos(ph_p) =',pharg pharg=1.0d0 else if (pharg.lt.-1.0d0) then write(*,*) ' phi_p error : cos(ph_p) =',pharg pharg=-1.0d0 end if php = dacos(pharg) sinphp = (dsin(thsc)*dcos(phsc)*dsin(ph0)*dcos(th0) + & dsin(thsc)*dsin(phsc)*dcos(ph0) + & dcos(thsc)*dsin(ph0)*dsin(th0) )/dsin(thp) if (sinphp.lt.0) php = 2.0d0*pi - php c temporarily not transforming thsc and phsc c if (j.eq.1) print *, " !Not transforming angles! (see target.f)" c thp = thsc c php = phsc c Cross section ......................... if(q.le.0.562803) then call intq(T, k0, q, thsc, xst) xstot = xst else c domga = phsc*(1.-dcos(thsc))/(2.*pi) ! don't need -- bec call sommers(T, en0, dxs) c xstot = domga*dxs/3.33 !! why /3.33?? bec xstot = dxs ! bec c We want probabilty of scattering, found from interaction length (dreact) found below. c This depends on the total cross section. The angular dependence is taken care of c in choosing q from S(q). To try and include some angular dependence here as well, c would be redundant. (see dsg0.f for same modification) -- bec endif c ............................................ c c the following lines are to turn off the cross section if you want to eliminate scattering effect c -- comment them out if you want scattering !! c if (j.eq.1) print *, " !Non-scattering run! (see target.f)" xstot = 1.d-6 c c call heden(T,hedens) denat = hedens*avog/he4mol call nden(T,norm) dreact = -1.0*log(vec(2))/(denat*xstot*norm) thxp = datan(dtan(thp)*dcos(php)) thyp = datan(dtan(thp)*dsin(php)) rfx = dreact*dsin(th0)*dcos(ph0) rfy = dreact*dsin(th0)*dsin(ph0) rfz = dreact*dcos(th0) rdist(1) = x0 + rfx rdist(2) = y0 + rfy rdist(3) = z0 + rfz scatflg=.false. if( (rdist(1).ge.tx1.and.rdist(1).le.tx2) .and. & (rdist(2).ge.ty1.and.rdist(2).le.ty2) .and. & (rdist(3).ge.z0.and.rdist(3).le.(z0+tz)) ) then !scat in tar xdis = rdist(1) + (tz-rfz)*dtan(thxp) ! x loc at end of tar ydis = rdist(2) + (tz-rfz)*dtan(thyp) ! y loc at end of tar scatflg = .true. if( (xdis.ge.tx1.and.xdis.le.tx2) .and. & (ydis.ge.ty1.and.ydis.le.ty2) ) then x = xdis y = ydis z = z0 + tz c thx = thxp c thy = thyp thx = 0.0d0 thy = 0.0d0 enp = ep else x = -1.d6 y = -1.d6 z = -1.d6 thx = -1.d6 thy = -1.d6 enp = -1.d6 endif else xdis = x0 + tz*dtan(thx0) ydis = y0 + tz*dtan(thy0) if( (xdis.ge.tx1.and.xdis.le.tx2) .and. & (ydis.ge.ty1.and.ydis.le.ty2) ) then x = xdis y = ydis z = z0 + tz thx = thx0 thy = thy0 enp = en0 else x = -1.d6 y = -1.d6 z = -1.d6 thx = -1.d6 thy = -1.d6 enp = -1.d6 endif endif c endif else x = -1.d6 y = -1.d6 z = -1.d6 thx = -1.d6 thy = -1.d6 enp = -1.d6 endif return end