program main_cwn implicit none include 'const.inc' integer i, j, nevts, counter, lim, Event, cin, cot, trig parameter(nevts=1e6) real xlam, HRNDM1, vect(5), lm(3), const, xinit, yinit, plr,xtest character*28 outfname double precision th0, ph0, thindex double precision en0, x0, y0, z0, thx0, thy0 double precision en1, x1, y1, z1, thx1, thy1 double precision en2, x2, y2, z2, thx2, thy2 double precision x3, y3, z3, thx3, thy3 double precision x4, y4, z4, thx4, thy4 double precision x5, y5, z5, thx5, thy5 double precision x6, y6, z6, thx6, thy6 double precision x7, y7, z7, thx7, thy7 double precision x8, y8, z8, thx8, thy8 double precision x9, y9, z9, thx9, thy9 double precision gl, gin, got integer ranno integer seedc integer time1, n1, n2, n3 C parameter(ranno=1) double precision T parameter (gl = 116.0, ! Guide length (cm) + gin = 0.001, ! Input critical index (rad/Ang) + T = 2.10) ! helium temp. (K) c + T = 4.25) ! helium temp. (K) c + got = 0.001) ! Output critical index (rad/Ang) double precision tx1, tx2, tx3, tx4, ty1, ty2, tz parameter (tx1 = -3.0, ! cm + tx2 = -0.35, ! cm + tx3 = 0.35, ! cm + tx4 = 3.0, ! cm + ty1 = -2.5, ! cm + ty2 = 2.5, ! cm + tz = 41.6) ! cm double precision gab1, gab2, gab3 parameter (gab1 = 9.0, ! gab between incoil and target (cm) + gab2 = 6.8, ! cm + gab3 = 12.0) ! cm double precision r2d, lambda double precision pl1_ns, th1_ns, pl2_ns, th2_ns, pl_ns, th_ns double precision scd1, scd2, pl1a_sc, pl1b_sc, th1a_sc, th1b_sc double precision pl2a_sc, pl2b_sc, th2a_sc, th2b_sc double precision pl1_sc, th1_sc, pl2_sc, th2_sc, pl3, th3 c----------------------------------- integer qmin, qmax double precision k0, qdummy, q c----------------------------------- 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) c time1=time() c ranno=mod(time1,2000000) c ranno=mod(time1,1089926518) print *,"asking for seed" read(*,*) ranno; c ranno=seedc() c print *, "ranNo is ",ranno, time1 c stop c stop/N/B/msarsour/spin/new_ih/ call hlimit(max_paw) iquest(10) = 6500000 c c----------------------------------------------------------------------- c----------------------------------------------------------------------- c call hrget(100,'lam.hbook','') c call hrget(102,'xyspot.hbook','') call hrget(151,'gradx.hbook','') call hrget(152,'gradx.hbook','') call hrget(301,'sq.hbook','') c----------------------------------------------------------------------- c read(*,*) ranno call rluxgo(3,ranno,0,0) ! RANLUX luxury level 2 cc call rluxgo(3,710,0,0) ! RANLUX luxury level 2 c----------------------------------------------------------------------- call hbook1(1001,'x (cm)',100,-3.5,3.3,0.) call hbook1(1002,'y (cm)',100,-3.,3.,0.) call hbook1(1003,'thx (mrad)',160,-1.6,1.6,0.) call hbook1(1004,'thy (mrad)',160,-1.6,1.6,0.) call hbook1(1005,'E (meV)/ non-scatt',100,0.,100.,0.) call hbook1(1006,'Path length (cm)/ non-scatt',100,41.5,41.7,0.) call hbook1(1007,'th_rotation (mrad)/ non-scatt',100,0.,1000.,0.) call hbook1(1008,'Path length (cm)/ non-scatt',100,41.5,41.7,0.) call hbook1(1009,'th_rotation (mrad)/ non-scatt',100,0.,1000.,0.) call hbook1(1010,'Path length (cm)/ non-scatt',100,41.5,41.7,0.) call hbook1(1011,'th_rotation (mrad)/ non-scatt',100,0.,1000.,0.) call hbook1(1012,'E (meV)/ scatt',100,0.,100.,0.) call hbook1(1013,'Path length (cm)/ scatt',100,41.,45.,0.) call hbook1(1014,'th_rotation (mrad)/ scatt',100,0.,1000.,0.) call hbook1(1015,'Path length (cm)/ scatt',100,41.,45.,0.) call hbook1(1016,'th_rotation (mrad)/ scatt',100,0.,1000.,0.) call hbook1(1017,'Path length (cm)/ scatt',100,41.,45.,0.) call hbook1(1018,'th_rotation (mrad)/ scatt',100,0.,1000.,0.) * call hbook1(2001,'x (cm)',100,-3.5,3.3,0.) call hbook1(2002,'y (cm)',100,-3.,3.,0.) call hbook1(2003,'thx (mrad)',160,-1.6,1.6,0.) call hbook1(2004,'thy (mrad)',160,-1.6,1.6,0.) call hbook1(2005,'E (meV)/ non-scatt',100,0.,100.,0.) call hbook1(2006,'Path length (cm)/ non-scatt',100,41.5,41.7,0.) call hbook1(2007,'th_rotation (mrad)/ non-scatt',100,0.,1000.,0.) call hbook1(2008,'Path length (cm)/ non-scatt',100,41.5,41.7,0.) call hbook1(2009,'th_rotation (mrad)/ non-scatt',100,0.,1000.,0.) call hbook1(2010,'Path length (cm)/ non-scatt',100,41.5,41.7,0.) call hbook1(2011,'th_rotation (mrad)/ non-scatt',100,0.,1000.,0.) call hbook1(2012,'E (meV)/ scatt',100,0.,100.,0.) call hbook1(2013,'Path length (cm)/ scatt',100,41.,45.,0.) call hbook1(2014,'th_rotation (mrad)/ scatt',100,0.,1000.,0.) call hbook1(2015,'Path length (cm)/ scatt',100,41.,45.,0.) call hbook1(2016,'th_rotation (mrad)/ scatt',100,0.,1000.,0.) call hbook1(2017,'Path length (cm)/ scatt',100,41.,45.,0.) call hbook1(2018,'th_rotation (mrad)/ scatt',100,0.,1000.,0.) * call hbook1(3001,'x (cm)',100,-3.5,3.3,0.) call hbook1(3002,'y (cm)',100,-3.,3.,0.) call hbook1(3003,'thx (mrad)',200,0.,0.1,0.) call hbook1(3004,'thy (mrad)',200,0.,0.1,0.) call hbook1(3005,'E (meV)',1000,0.,100.,0.) call hbook1(4005,'Lambda (Angs)',1000,0.,100.,0.) call hbook1(3006,'Path length (cm)',200,83.15,83.25,0.) call hbook1(3007,'th_rotation (mrad)',200,0.,1000.,0.) call hbook1(3008,'Path length (cm)',200,83.15,83.25,0.) call hbook1(3009,'th_rotation (mrad)',200,0.,1000.,0.) call hbook1(3010,'Path length (cm)',200,83.15,83.25,0.) call hbook1(3011,'th_rotation (mrad)',200,0.,1000.,0.) call hbook1(3012,'E (meV)/ non-scatt',200,0.,100.,0.) call hbook1(3013,'E (meV)/ non-scatt',200,0.,100.,0.) call hbook1(3014,'E (meV)/ non-scatt',200,0.,100.,0.) c c------------Input file ------------------------- c open(39,file='inp.dat',status='old') c read(39,*) gl, gin, T c------------------------------------------------ r2d=180./pi c------------ counter = 0 Event = 0 cin=0 cot=0 c got = gin c do 100 i = 1, nevts c Event = 1 c c Generate initial neutron wavelength c xlam = HRNDM1(100) c c The energy calculated from the neutron's wavelength c if(xlam.lt.1.) goto 100 en0 = (1/(2*m_n))*(2*pi*h_par/xlam)**2 c c Generates the initial trajectory of the neutron. c call ranlux(vect,4) thindex = 0.003*xlam th0 = thindex*vect(1) ! polar angle ph0 = 2.00*pi*vect(2) ! azimuthal angle (radians) c call HRNDM2(102,xinit,yinit) c x0=xinit/10. ! cm c y0=yinit/10. ! cm c write(18,*) x0, y0; goto 100 xtest = 2.75*(2.*vect(3)-1.) if(xtest.gt.0.) then x0 = HRNDM1(151) elseif(xtest.lt.0) then x0 = HRNDM1(151) - 3. endif y0 = 2.25*(2.*vect(4)-1.) z0= 0.0d0 c c Calculate the incident angles in the xz-plane and yz-plane c thx0 = datan(dtan(th0)*dcos(ph0)) thy0 = datan(dtan(th0)*dsin(ph0)) c k0 = sqrt(2*m_n*en0)/h_par qmin = 0 qmax = int(k0*25.) ! for sq_1.hbook if(k0.gt.3.6) qmax=int(3.6*25.) if(i.gt.1) call hdelet(302) call hcopyr(301,302,"",qmin,qmax,0.,0.,"") qdummy=HRNDM1(301) q=HRNDM1(302) ******************************************************************** ******************************************************************** * Neutrons passing thruogh the guide * call guide(gl,gin,en0,x0,y0,z0,thx0,thy0,x1,y1,z1,thx1,thy1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 * ******************************************************************** ******************************************************************** * Neutrons passing thruogh the input coil * call incoil(gin,en0,x1,y1,z1,thx1,thy1,x2,y2,z2,thx2,thy2) if(x2.lt.-1d5.or.y2.lt.-1d5) goto 100 * ******************************************************************** ******************************************************************** * Neutrons passing thruogh the gab (9.0 cm) between the input * coil and the target. * x3 = x2 + gab1*dtan(thx2) y3 = y2 + gab1*dtan(thy2) z3 = z2 + gab1 thx3 = thx2 thy3 = thy2 * ******************************************************************** ******************************************************************** * Neutrons passing thruogh the 1st half of the target (41.6 cm). * if( (x3.ge.tx1.and.x3.le.tx2) ) then call empty(en0, x3, y3, z3, thx3, thy3, tx1, tx2, & ty1, ty2, tz, en1, x4, y4, z4, thx4, thy4) elseif( (x3.ge.tx3.and.x3.le.tx4) ) then call target(q,T, en0, x3, y3, z3, thx3, thy3, tx3, tx4, & ty1, ty2, tz, en1, x4, y4, z4, thx4, thy4) else x4 = -1.d6 y4 = -1.d6 z4 = -1.d6 thx4 = -1.d6 thy4 = -1.d6 en1 = -1.d6 endif * if(x4.lt.-1d5.or.y4.lt.-1d5) goto 100 * call hf1(1001,real(x4),1.) call hf1(1002,real(y4),1.) call hf1(1003,real(thx4),1.) call hf1(1004,real(thy4),1.) if(en0.eq.en1) then pl1_ns = sqrt((x4-x3)**2+(y4-y3)**2+(z4-z3)**2) th1_ns = gy*mag_B*pl1_ns*0.1/sqrt(2.*en1/m_n) c write(16,*) th1_ns pl1_sc = 0. th1_sc = 0. call hf1(1005,real(en1),1.) call hf1(1006,real(pl1_ns),1.) call hf1(1007,real(th1_ns),1.) if(x3.gt.0..and.x4.gt.0.) then call hf1(1008,real(pl1_ns),1.) call hf1(1009,real(th1_ns),1.) elseif(x3.lt.0..and.x4.lt.0.) then call hf1(1010,real(pl1_ns),1.) call hf1(1011,real(th1_ns),1.) endif else scd1 = (x4-x3-(tz*dtan(thx4)))/(dtan(thx3)-dtan(thx4)) pl1a_sc = sqrt((scd1*dtan(thx3))**2+(scd1*dtan(thy3))**2+ & scd1**2) pl1b_sc = sqrt( ((tz-scd1)*dtan(thx4))**2 & + ((tz-scd1)*dtan(thy4))**2+(tz-scd1)**2) th1a_sc = gy*mag_B*pl1a_sc*0.1/sqrt(2.*en0/m_n) th1b_sc = gy*mag_B*pl1b_sc*0.1/sqrt(2.*en1/m_n) c write(17,*) th1a_sc,th1b_sc pl1_ns = 0. th1_ns = 0. pl1_sc = pl1a_sc + pl1b_sc th1_sc = th1a_sc + th1b_sc call hf1(1012,real(en1),1.) call hf1(1013,real(pl1_sc),1.) call hf1(1014,real(th1_sc),1.) if(x3.gt.0..and.x4.gt.0.) then call hf1(1015,real(pl1_sc),1.) call hf1(1016,real(th1_sc),1.) elseif(x3.lt.0..and.x4.lt.0.) then call hf1(1017,real(pl1_sc),1.) call hf1(1018,real(th1_sc),1.) endif endif * ******************************************************************** ******************************************************************** * Neutrons passing thruogh the gab (pi-coil = 7.3 cm) between the * two tagets. * x5 = x4 + gab2*dtan(thx4) y5 = y4 + gab2*dtan(thy4) z5 = z4 + gab2 thx5 = thx4 thy5 = thy4 * ******************************************************************** ******************************************************************** * Neutrons passing thruogh the 2nd part of the target (41.6 cm). * c call hf1(203,real(x5),1.) if( (x5.gt.tx1.and.x5.lt.tx2) ) then call target(q,T, en1, x5, y5, z5, thx5, thy5, tx1, tx2, & ty1, ty2, tz, en2, x6, y6, z6, thx6, thy6) elseif( (x5.gt.tx3.and.x5.lt.tx4) ) then call empty(en1, x5, y5, z5, thx5, thy5, tx3, tx4, & ty1, ty2, tz, en2, x6, y6, z6, thx6, thy6) else x6 = -1.d6 y6 = -1.d6 z6 = -1.d6 endif * if(x6.lt.-1d5.or.y6.lt.-1d5) goto 100 * call hf1(2001,real(x6),1.) call hf1(2002,real(y6),1.) call hf1(2003,real(thx6),1.) call hf1(2004,real(thy6),1.) if(en2.eq.en1) then pl2_ns = sqrt((x6-x5)**2+(y6-y5)**2+(z6-z5)**2) th2_ns = gy*mag_B*pl2_ns*0.1/sqrt(2.*en2/m_n) c write(18,*) th2_ns pl2_sc = 0. th2_sc = 0. call hf1(2005,real(en2),1.) call hf1(2006,real(pl2_ns),1.) call hf1(2007,real(th2_ns),1.) if(x5.gt.0..and.x6.gt.0.) then call hf1(2008,real(pl2_ns),1.) call hf1(2009,real(th2_ns),1.) elseif(x5.lt.0..and.x6.lt.0.) then call hf1(2010,real(pl2_ns),1.) call hf1(2011,real(th2_ns),1.) endif else scd2 = (x6-x5-(tz*dtan(thx6)))/(dtan(thx5)-dtan(thx6)) pl2a_sc = sqrt((scd2*dtan(thx5))**2+(scd2*dtan(thy5))**2+ & scd2**2) pl2b_sc = sqrt( ((tz-scd2)*dtan(thx6))**2 & + ((tz-scd2)*dtan(thy6))**2+(tz-scd2)**2) th2a_sc = gy*mag_B*pl2a_sc*0.1/sqrt(2.*en1/m_n) th2b_sc = gy*mag_B*pl2b_sc*0.1/sqrt(2.*en2/m_n) write(19,*) th2a_sc,th2b_sc pl2_ns = 0. th2_ns = 0. pl2_sc = pl2a_sc + pl2b_sc th2_sc = th2a_sc + th2b_sc call hf1(2012,real(en2),1.) call hf1(2013,real(pl2_sc),1.) call hf1(2014,real(th2_sc),1.) if(x5.gt.0..and.x6.gt.0.) then call hf1(2015,real(pl2_sc),1.) call hf1(2016,real(th2_sc),1.) elseif(x5.lt.0..and.x6.lt.0.) then call hf1(2017,real(pl2_sc),1.) call hf1(2018,real(th2_sc),1.) endif endif * ******************************************************************** ******************************************************************** * Neutrons passing thruogh the gab (9.6 cm) between the taget and * the output coil. * x7 = x6 + gab3*dtan(thx6) y7 = y6 + gab3*dtan(thy6) z7 = z6 + gab3 thx7 = thx6 thy7 = thy6 * ******************************************************************** ******************************************************************** * Neutrons passing thruogh the output coil. * call outcoil(got,en2,x7,y7,z7,thx7,thy7,x8,y8,z8,thx8,thy8) if(x8.lt.-1d5.or.y8.lt.-1d5) goto 100 * ******************************************************************** ******************************************************************** * Neutrons passing thruogh the supper mirror. * call sm(en2,x8,y8,z8,thx8,thy8,x9,y9,z9,thx9,thy9) if(x9.lt.-1d5.or.y9.lt.-1d5) goto 100 cot=cot+1 * ******************************************************************** ******************************************************************** * call hf1(3001,real(x9),1.) call hf1(3002,real(y9),1.) call hf1(3003,real(thx9),1.) call hf1(3004,real(thy9),1.) pl3 = pl1_ns + pl2_ns + pl1_sc + pl2_sc th3 = th1_ns + th2_ns + th1_sc + th2_sc c write(20,*) th3 if(en0.ne.en1 .or. en1.ne.en2) & write(*,*) pl1_ns, pl2_ns, pl1_sc, pl2_sc call hf1(3005,real(en2),1.) lambda=(2.*pi*h_par)/sqrt(2*m_n*en2) call hf1(4005,real(lambda),1.) call hf1(3006,real(pl3),1.) call hf1(3007,real(th3),1.) if(x3.gt.0..and.x6.gt.0.) then call hf1(3008,real(pl3),1.) call hf1(3009,real(th3),1.) call hf1(3012,real(en2),1.) elseif(x3.lt.0..and.x6.lt.0.) then call hf1(3010,real(pl3),1.) call hf1(3011,real(th3),1.) call hf1(3013,real(en2),1.) endif if(en0.ne.en1 .or. en1.ne.en2) then call hf1(3014,real(en2),1.) endif ******************************************************************** * 34 format(i12,2x,a6) do j=1,100 lim = j*1e7 if(i.eq.lim) write(*,34) i, 'events' enddo 100 continue write(*,*) cin,' ', cot c c ------------------------- if(ranno<1e1) then Write(outfname,'(1Hh,I1,6H.hbook)')ranno elseif(ranno<1e2) then Write(outfname,'(1Hh,I2,6H.hbook)')ranno elseif(ranno<1e3) then Write(outfname,'(1Hh,I3,6H.hbook)')ranno elseif(ranno<1e4) then Write(outfname,'(1Hh,I4,6H.hbook)')ranno elseif(ranno<1e5) then Write(outfname,'(1Hh,I5,6H.hbook)')ranno elseif(ranno<1e6) then Write(outfname,'(1Hh,I6,6H.hbook)')ranno elseif(ranno<1e7) then Write(outfname,'(1Hh,I7,6H.hbook)')ranno elseif(ranno<1e8) then Write(outfname,'(1Hh,I8,6H.hbook)')ranno endif call hrput(0,outfname,'N') c ------------------------- c return end