program main implicit none include 'const.inc' integer i, j, n, nevts, counter, ranno parameter(nevts=10) parameter(n=19) real lm, HRNDM1, vect(6), const, xinit, yinit, plr,xtest double precision th0, ph0, dz, x(n), y(n), z(n) double precision bx(n), by(n), bz(n), px(n), py(n), pz(n) double precision thxy(n),thxz(n),thyx(n),thyz(n),thzx(n),thzy(n) double precision thx(n), thy(n), thz(n), phx(n), phy(n), phz(n) double precision lmin, lmax, thindex double precision xmin, xmax, ymin, ymax, zmin, zmax double precision throtx, throty, throtz, en double precision xd1, xd2, yd1, yd2, t1, t2, dt parameter(xd1=-3.5, xd2=3.5) parameter(yd1=-2.5, yd2=2.5) c c ------------------------- c fill the ntuple 'spin.rz' //// ntuple variables .... integer max_paw, iquest integer istat, icycle parameter(max_paw = 10000000) ! PAW common block size real quest, p common/quest/iquest(100) common/pawc/p(max_paw) real lam, xi, yi, zi, thi, phi, xf, yf, zf real polx, poly, polz, rotxy, rotxz, rotyx, rotyz, rotzx, rotzy common /spinblk/ i,lam, & xi, yi, zi, thi, phi, & xf, yf, zf, polx, poly, polz, & rotxy, rotxz, rotyx, rotyz, rotzx, rotzy call hlimit(max_paw) iquest(10) = 6500000 call hropen(1,'Myntuple','ntbx0by0.rz','NQ',8190,istat) call hbnt(10,'spinnt',' ') c-------------------------------- call hbname(10,'spin',i,'i:I,lam:R,xi:R,yi:R,zi:R,'// & 'thi:R, phi:R, xf:R, yf:R, zf:R, polx:R, poly:R, polz:R,'// & 'rotxy:R, rotxz:R, rotyx:R, rotyz:R, rotzx:R, rotzy:R') c----------------------------------------------------------------------- call hrget(100,'lam.hbook','') c----------------------------------------------------------------------- c t1=time() ccc ranno=mod(t1,1576879653) ccc print *,"asking for seed" ccc read(*,*) ranno; c c call rluxgo(3,ranno,0,0) ! RANLUX luxury level 3 call rluxgo(3,11,0,0) ! RANLUX luxury level 3 c----------------------------------------------------------------------- c lmin = 2.0 ; lmax = 15.0 xmin = -3.5 ; xmax = 3.5 ymin = -2.0 ; ymax = 2.0 zmin = -9.0 ; zmax = 9.0 c do 100 i = 1, nevts c lm = HRNDM1(100) c lm=5 c call ranlux(vect,4) x(1) = xmin + (xmax-xmin)*vect(1) y(1) = ymin + (ymax-ymin)*vect(2) z(1) = zmin thindex = 0.003*lm th0 = thindex*vect(3) ph0 = 2.00*pi*vect(4) px(1) = 0.0 py(1) = 1.0 pz(1) = 0.0 dz=(zmax-zmin+1)/n const=gy*m_n*lm*dz*1d8/h_blanck counter=1 do 200 j=2,n counter=counter+1 z(j) = zmin+(j-1)*dz x(j) = x(j-1)+z(j)*dsin(th0)*dcos(ph0) y(j) = y(j-1)+z(j)*dsin(th0)*dsin(ph0) c write(*,*) i,j,x(j),y(j),z(j); goto 200 if(x(j).lt.xd1.or.x(j).gt.xd2) goto 100 if(y(j).lt.yd1.or.y(j).gt.yd2) goto 100 c call mag_b(x(j),y(j),z(j),bx(j),48,1) call mag_b(x(j),y(j),z(j),by(j),48,2) call mag_b(x(j),y(j),z(j),bz(j),60,3) c thxy(j)=const*(px(j-1)*bz(j)) thxz(j)=const*(px(j-1)*by(j)) phx(j)=datan(dtan(thxz(j))/dtan(thxy(j))) thx(j)=datan(dtan(thxy(j))/dcos(phx(j))) if(thxy(j).eq.0..and.thxy(j).eq.0.) then thx(j)=0. phx(j)=0. elseif(thxy(j).eq.0.) then thx(j)=thxz(j) phx(j)=pi/2. elseif(thxz(j).eq.0.) then thx(j)=thxy(j) phx(j)=0. endif c thyx(j)=const*(py(j-1)*bz(j)) thyz(j)=const*(py(j-1)*bx(j)) phy(j)=datan(dtan(thyx(j))/dtan(thyz(j))) thy(j)=datan(dtan(thyz(j))/dcos(phy(j))) if(thyx(j).eq.0..and.thyz(j).eq.0.) then thy(j)=0. phy(j)=0. elseif(thyz(j).eq.0.) then thy(j)=thyx(j) phy(j)=pi/2. elseif(thyx(j).eq.0.) then thy(j)=thyz(j) phy(j)=0. endif c thzx(j)=const*(pz(j-1)*by(j)) thzy(j)=const*(pz(j-1)*bx(j)) phz(j)=datan(dtan(thzy(j))/dtan(thzx(j))) thz(j)=datan(dtan(thzx(j))/dcos(phz(j))) if(thzx(j).eq.0..and.thzy(j).eq.0.) then thz(j)=0. phz(j)=0. elseif(thzx(j).eq.0.) then thz(j)=thzy(j) phz(j)=pi/2. elseif(thzy(j).eq.0.) then thz(j)=thzx(j) phz(j)=0. endif c px(j) = px(j-1)*dcos(thx(j)) + & py(j-1)*dsin(thy(j))*dsin(phy(j)) + & pz(j-1)*dsin(thz(j))*dcos(phz(j)) py(j) = px(j-1)*dsin(thx(j))*dcos(phx(j)) + & py(j-1)*dcos(thy(j)) + & pz(j-1)*dsin(thz(j))*dsin(phz(j)) pz(j) = px(j-1)*dsin(thx(j))*dsin(phx(j)) + & py(j-1)*dsin(thy(j))*dcos(phy(j)) + & pz(j-1)*dcos(thz(j)) c 200 continue c c c Fill the ntuple variables c lam=lm xi = x(1) yi = y(1) zi = z(1) thi = th0 phi = ph0 xf = x(counter) yf = y(counter) zf = z(counter) polx = px(counter) poly = py(counter) polz = pz(counter) rotxy = datan(py(counter)/px(counter)) rotxz = datan(pz(counter)/px(counter)) rotyx = datan(px(counter)/py(counter)) rotyz = datan(pz(counter)/py(counter)) rotzx = datan(px(counter)/pz(counter)) rotzy = datan(py(counter)/pz(counter)) c write(*,*) polx, poly, polz call hfnt(10) c write(*,*) counter, z(counter),zi,zf c stop 100 continue c t2=time() dt=t2-t1 write(*,*) dt," seconds" c c ------------------------- call hrout(10,icycle,' ') call hrend('Myntuple') c ------------------------- c return end