42c42,43 < double precision en1, x1, y1, z1, thx1, thy1 --- > double precision en1, x1, y1, z1, thx1, thy1 > real xtmp, ytmp 198c199,200 < cbec10 call hrget(700,'hbooks/sqw.hbook','') ! New S(q,w) --- > cbec10 call hrget(700,'hbooks/sqw.hbook','') ! New S(q,w) > call hrget(200,'hbooks/psm.hbook','') !x-y out of PSM -bec 226,235c228,244 < call hbook2(1110,'x-y PSM out',55,-5.5,5.5,55, < & -5.5,5.5,0.) !-bec < call hbook2(1111,'x-y incoil out',55,-5.5,5.5,55, < & -5.5,5.5,0.) !-bec < call hbook2(2222,'x-y target out',55,-5.5,5.5,55, < & -5.5,5.5,0.) !-bec < call hbook2(3333,'x-y outcoil out',55,-5.5,5.5,55, < & -5.5,5.5,0.) !-bec < call hbook2(4444,'x-y Detector',55,-5.5,5.5,55, < & -5.5,5.5,0.) !-bec --- > call hbook2(2220,'Flux(xy) PSM out',300,-2.98,3.0,300, > & -2.98,3.0,0.) !-bec > call hbook2(2221,'Flux(xy) incoil out',300,-2.98,3.0,300, > & -2.98,3.0,0.) !-bec > call hbook2(2222,'Flux(xy) target out',350,-3.48,3.5,350, > & -3.48,3.5,0.) !-bec > call hbook2(3333,'Flux(xy) outcoil out',350,-3.48,3.5,350, > & -3.48,3.5,0.) !-bec > call hbook2(4444,'Flux(xy) Detector',350,-3.48,3.5,350, > & -3.48,3.5,0.) !-bec > call hbook1(441,'Flux(x) PSM (/s)',300,-3.0,3.0,0.) > call hbook1(442,'Flux(y) PSM (/s)',300,-3.0,3.0,0.) > call hbook1(443,'Flux(x) InGuide (/s)',300,-3.0,3.0,0.) > call hbook1(444,'Flux(y) InGuide (/s)',300,-3.0,3.0,0.) > call hbook1(445,'Flux(x) OutGuide (/s)',350,-3.5,3.5,0.) > call hbook1(446,'Flux(y) OutGuide (/s)',350,-3.5,3.5,0.) > 327,330c336,339 < call hbook1(6005,'Path length (cm) Det scat West',300,492.20, < & 492.60,0.) < call hbook1(6007,'Path length (cm) Det scat East',300,492.20, < & 492.60,0.) --- > call hbook1(6005,'Path length (cm) Det scat West',300,492.50, > & 493.50,0.) > call hbook1(6007,'Path length (cm) Det scat East',300,492.50, > & 493.50,0.) 376,396c385,405 < call hbook2(6030,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40., < & 300,-800.,0.,0.) < call hbook2(6031,'th_rotsc (mrad) vs. E(meV) det East',800,0.,40., < & 300,-800.,0.,0.) < call hbook2(6020,'th_rot (mrad) vs. E (meV) det West',800,0.,40., < & 300,-800.,0.,0.) < call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300, < & 480.50,480.90,300,-800.,0.,0.) < call hbook2(6022,'th_rot (mrad) vs. E (meV) det East',800,0.,40., < & 300,-800.,0.,0.) < call hbook2(6023,'th rot (mrad) vs. Path len (cm) det East',300, < & 480.50,480.90,300,-800.,0.,0.) < call hbook2(6040,'Cr-ov. th_rot (mrad) vs. E (meV) det West', < & 800,0.,40.,300,-800.,0.,0.) < call hbook2(6041,'Cr-ov. th rot (mrad) vs. Path len (cm) det West', < & 300,480.50,480.90,300,-800.,0.,0.) < call hbook2(6042,'Cr-ov. th_rot (mrad) vs. E (meV) det East', < & 800,0.,40.,300,-800.,0.,0.) < call hbook2(6043,'Cr-ov. th rot (mrad) vs. Path len (cm) det East', < & 300,480.50,480.90,300,-800.,0,0.) < --- > c call hbook2(6030,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40., > c & 300,-800.,0.,0.) > c call hbook2(6031,'th_rotsc (mrad) vs. E(meV) det East',800,0.,40., > c & 300,-800.,0.,0.) > c call hbook2(6020,'th_rot (mrad) vs. E (meV) det West',800,0.,40., > c & 300,-800.,0.,0.) > c call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300, > c & 480.50,480.90,300,-800.,0.,0.) > c call hbook2(6022,'th_rot (mrad) vs. E (meV) det East',800,0.,40., > c & 300,-800.,0.,0.) > c call hbook2(6023,'th rot (mrad) vs. Path len (cm) det East',300, > c & 480.50,480.90,300,-800.,0.,0.) > c call hbook2(6040,'Cr-ov. th_rot (mrad) vs. E (meV) det West', > c & 800,0.,40.,300,-800.,0.,0.) > c call hbook2(6041,'Cr-ov. th rot (mrad) vs. Path len (cm) det West', > c & 300,480.50,480.90,300,-800.,0.,0.) > c call hbook2(6042,'Cr-ov. th_rot (mrad) vs. E (meV) det East', > c & 800,0.,40.,300,-800.,0.,0.) > c call hbook2(6043,'Cr-ov. th rot (mrad) vs. Path len (cm) det East', > c & 300,480.50,480.90,300,-800.,0,0.) > 420,440c429,449 < call hbook2(6030,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40., < & 300,-800.,0.,0.) < call hbook2(6031,'th_rotsc (mrad) vs. E(meV) det East',800,0.,40., < & 300,-800.,0.,0.) < call hbook2(6020,'th_rot (mrad) vs. E (meV) det West',800,0.,40., < & 300,0.,60.,0.) < call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300, < & 480.50,480.90,300,0.,60.,0.) < call hbook2(6022,'th_rot (mrad) vs. E (meV) det East',800,0.,40., < & 300,0.,60.,0.) < call hbook2(6023,'th rot (mrad) vs. Path len (cm) det East',300, < & 480.50,480.90,300,0.,60.,0.) < call hbook2(6040,'Cr-ov. th_rot (mrad) vs. E (meV) det West', < & 800,0.,40.,300,0.,60.,0.) < call hbook2(6041,'Cr-ov. th rot (mrad) vs. Path len (cm) det West', < & 300,480.50,480.90,300,0.,60.,0.) < call hbook2(6042,'Cr-ov. th_rot (mrad) vs. E (meV) det East', < & 800,0.,40.,300,0.,60.,0.) < call hbook2(6043,'Cr-ov. th rot (mrad) vs. Path len (cm) det East', < & 300,480.50,480.90,300,0.,60,0.) < --- > c call hbook2(6030,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40., > c & 300,-800.,0.,0.) > c call hbook2(6031,'th_rotsc (mrad) vs. E(meV) det East',800,0.,40., > c & 300,-800.,0.,0.) > c call hbook2(6020,'th_rot (mrad) vs. E (meV) det West',800,0.,40., > c & 300,0.,60.,0.) > c call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300, > c & 480.50,480.90,300,0.,60.,0.) > c call hbook2(6022,'th_rot (mrad) vs. E (meV) det East',800,0.,40., > c & 300,0.,60.,0.) > c call hbook2(6023,'th rot (mrad) vs. Path len (cm) det East',300, > c & 480.50,480.90,300,0.,60.,0.) > c call hbook2(6040,'Cr-ov. th_rot (mrad) vs. E (meV) det West', > c & 800,0.,40.,300,0.,60.,0.) > c call hbook2(6041,'Cr-ov. th rot (mrad) vs. Path len (cm) det West', > c & 300,480.50,480.90,300,0.,60.,0.) > c call hbook2(6042,'Cr-ov. th_rot (mrad) vs. E (meV) det East', > c & 800,0.,40.,300,0.,60.,0.) > c call hbook2(6043,'Cr-ov. th rot (mrad) vs. Path len (cm) det East', > c & 300,480.50,480.90,300,0.,60,0.) > 468,488c477,497 < call hbook2(6030,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40., < & 300,-5.,50.,0.) < call hbook2(6031,'th_rotsc (mrad) vs. E(meV) det East',800,0.,40., < & 300,-5.,50.,0.) < call hbook2(6020,'th_rot (mrad) vs. E (meV) det West',800,0.,40., < & 300,-10.,10.,0.) < call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300, < & 480.50,480.90,300,-10.,10.,0.) < call hbook2(6022,'th_rot (mrad) vs. E (meV) det East',800,0.,40., < & 300,-10.,10.,0.) < call hbook2(6023,'th rot (mrad) vs. Path len (cm) det East',300, < & 480.50,480.90,300,-10.,10.,0.) < call hbook2(6040,'Cr-ov. th_rot (mrad) vs. E (meV) det West', < & 800,0.,40.,300,-10.,10.,0.) < call hbook2(6041,'Cr-ov. th rot (mrad) vs. Path len (cm) det West', < & 300,480.50,480.90,300,-10.,10.,0.) < call hbook2(6042,'Cr-ov. th_rot (mrad) vs. E (meV) det East', < & 800,0.,40.,300,-10.,10.,0.) < call hbook2(6043,'Cr-ov. th rot (mrad) vs. Path len (cm) det East', < & 300,480.50,480.90,300,-10.,10,0.) < --- > c call hbook2(6030,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40., > c & 300,-5.,50.,0.) > c call hbook2(6031,'th_rotsc (mrad) vs. E(meV) det East',800,0.,40., > c & 300,-5.,50.,0.) > c call hbook2(6020,'th_rot (mrad) vs. E (meV) det West',800,0.,40., > c & 300,-10.,10.,0.) > c call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300, > c & 480.50,480.90,300,-10.,10.,0.) > c call hbook2(6022,'th_rot (mrad) vs. E (meV) det East',800,0.,40., > c & 300,-10.,10.,0.) > c call hbook2(6023,'th rot (mrad) vs. Path len (cm) det East',300, > c & 480.50,480.90,300,-10.,10.,0.) > c call hbook2(6040,'Cr-ov. th_rot (mrad) vs. E (meV) det West', > c & 800,0.,40.,300,-10.,10.,0.) > c call hbook2(6041,'Cr-ov. th rot (mrad) vs. Path len (cm) det West', > c & 300,480.50,480.90,300,-10.,10.,0.) > c call hbook2(6042,'Cr-ov. th_rot (mrad) vs. E (meV) det East', > c & 800,0.,40.,300,-10.,10.,0.) > c call hbook2(6043,'Cr-ov. th rot (mrad) vs. Path len (cm) det East', > c & 300,480.50,480.90,300,-10.,10,0.) > 513,533c522,542 < call hbook2(6030,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40., < & 300,-5.,50.,0.) < call hbook2(6031,'th_rotsc (mrad) vs. E(meV) det East',800,0.,40., < & 300,-5.,50.,0.) < call hbook2(6020,'th_rot (mrad) vs. E (meV) det West',800,0.,40., < & 300,-5.,50.,0.) < call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300, < & 480.50,480.90,300,-5.,50.,0.) < call hbook2(6022,'th_rot (mrad) vs. E (meV) det East',800,0.,40., < & 300,-5.,50.,0.) < call hbook2(6023,'th rot (mrad) vs. Path len (cm) det East',300, < & 480.50,480.90,300,-5.,50.,0.) < call hbook2(6040,'Cr-ov. th_rot (mrad) vs. E (meV) det West', < & 800,0.,40.,300,-5.,50.,0.) < call hbook2(6041,'Cr-ov. th rot (mrad) vs. Path len (cm) det West', < & 300,480.50,480.90,300,-5.,50.,0.) < call hbook2(6042,'Cr-ov. th_rot (mrad) vs. E (meV) det East', < & 800,0.,40.,300,-5.,50.,0.) < call hbook2(6043,'Cr-ov. th rot (mrad) vs. Path len (cm) det East', < & 300,480.50,480.90,300,-5.,50,0.) < --- > c call hbook2(6030,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40., > c & 300,-5.,50.,0.) > c call hbook2(6031,'th_rotsc (mrad) vs. E(meV) det East',800,0.,40., > c & 300,-5.,50.,0.) > c call hbook2(6020,'th_rot (mrad) vs. E (meV) det West',800,0.,40., > c & 300,-5.,50.,0.) > c call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300, > c & 480.50,480.90,300,-5.,50.,0.) > c call hbook2(6022,'th_rot (mrad) vs. E (meV) det East',800,0.,40., > c & 300,-5.,50.,0.) > c call hbook2(6023,'th rot (mrad) vs. Path len (cm) det East',300, > c & 480.50,480.90,300,-5.,50.,0.) > c call hbook2(6040,'Cr-ov. th_rot (mrad) vs. E (meV) det West', > c & 800,0.,40.,300,-5.,50.,0.) > c call hbook2(6041,'Cr-ov. th rot (mrad) vs. Path len (cm) det West', > c & 300,480.50,480.90,300,-5.,50.,0.) > c call hbook2(6042,'Cr-ov. th_rot (mrad) vs. E (meV) det East', > c & 800,0.,40.,300,-5.,50.,0.) > c call hbook2(6043,'Cr-ov. th rot (mrad) vs. Path len (cm) det East', > c & 300,480.50,480.90,300,-5.,50,0.) > 687,690c696 < xtest = xin*(2.*vect(3)-1.) !-KG < if (i.eq.1) then < write(*,*) ' ! Running with No Position (x) gradient !' < write(3,*) ' ! Running with No Position (x) gradient !' --- > if (i.eq.1) then 695,721d700 < end if < < x0 = xtest < xstart = x0 !-KG < y0 = yin*(2.*vect(4)-1.) < < if(x0.gt.xmax) xmax = x0 < if(y0.gt.ymax) ymax = y0 < < ccccccccc gradient < cc < cc x0 = 3.0d0 - 2.d0*HRNDM1(151) !bec R to L gradient < cc < cc x0 = 2.d0*HRNDM1(151) - 3.0d0 !bec L to R gradient < cc < cc this gives double gradient -- gradient across each half... < cc if(xtest.gt.0.) then < cc x0 = HRNDM1(151) ! gives gradient only acrros half < ccc x0 = 3.d0 - HRNDM1(151) < cc elseif(xtest.lt.0) then < cc x0 = HRNDM1(151) - 3. < ccc x0 = -HRNDM1(151) < cc endif < cc < cc xstart = x0 < ccccccccc < 723d701 < if (i.eq.1) then 748,754d725 < if (sptmflg) then < write(6,*) ' ! Septum Upstream ! ' < write(3,*) ' ! Septum Upstream ! ' < else < write(6,*) ' ! No Septum Upstream ! ' < write(3,*) ' ! No Septum Upstream ! ' < end if 758a730,774 > ccccccccccccccccccccccc No x-y gradient > c if (i.eq.1) then > c write(*,*) ' ! Running with No Position (x) gradient !' > c write(3,*) ' ! Running with No Position (x) gradient !' > c end if > c xtest = xin*(2.*vect(3)-1.) !-KG > c x0 = xtest > c xstart = x0 !-KG > c y0 = yin*(2.*vect(4)-1.) > c > ccccccccc gradient > cc > cc x0 = 3.0d0 - 2.d0*HRNDM1(151) !bec R to L gradient > cc > cc x0 = 2.d0*HRNDM1(151) - 3.0d0 !bec L to R gradient > cc > cc this gives double gradient -- gradient across each half... > cc if(xtest.gt.0.) then > cc x0 = HRNDM1(151) ! gives gradient only acrros half > ccc x0 = 3.d0 - HRNDM1(151) > cc elseif(xtest.lt.0) then > cc x0 = HRNDM1(151) - 3. > ccc x0 = -HRNDM1(151) > cc endif > cc > cc xstart = x0 > ccccccccc > > cccccccccccc x-y from PSM image > c > call HRNDM2(200,xtmp,ytmp) > x0= dble(xtmp) > y0 = dble(ytmp) > xtest = x0 > xstart = x0 > > > if ((abs(x0).gt.xin).or.(abs(y0).gt.yin)) goto 100 > > if (x0**2+y0**2.le.0.625**2) > & INPKflux = INPKflux + 1/pi/0.625**2/nevts > > > > cccccccccccc 762,772c778,790 < c if (i.eq.1) then < c write(6,*) ' !!! Inverted Parabola x-y dist!!' < c write(3,*) ' !!! Inverted Parabola x-y dist!!' < c end if < c < c xtest = HRNDM1(161) < c y0 = HRNDM1(162) < c x0 = xtest < c xstart = xtest < c if (x0.gt.xmax) xmax = x0 < c if (y0.gt.ymax) ymax = y0 --- > c if (i.eq.1) then > c write(6,*) ' !!! Inverted Parabola x-y dist!!' > c write(3,*) ' !!! Inverted Parabola x-y dist!!' > c end if > c > c xtest = HRNDM1(161) > c y0 = HRNDM1(162) > c x0 = xtest > c xstart = xtest > c if (x0.gt.xmax) xmax = x0 > c if (y0.gt.ymax) ymax = y0 > c > c if ((abs(x0).gt.xin).or.(abs(y0).gt.yin)) goto 100 774c792 < c if (x0**2+y0**2.le.0.625**2) --- > c if (x0**2+y0**2.le.0.625**2) 776,778c794,795 < c < c < c if ((abs(x0).gt.xin).or.(abs(y0).gt.yin)) goto 100 --- > c > c 797,809c814,826 < c thindex = ng6indexx*xlam ! NG6 about 1.7mrad/ang < c ph0 = 2.0d0*pi*vect(2) ! azimuthal angle (radians) < c < c if (thindex.lt.1E-7) thindex = 1E-7 < c if (abs(1.d0-(1.d0-dcos(thindex))*vect(1)).gt.1.0d0) then < c write(6,*) ' Arccos(theta_start)>1 ' < c goto 321 < c end if < c < c th0 = dacos(1.d0-(1.d0-dcos(thindex))*vect(1)) ! isotropic (0,thindex) bec < c thx0 = ratan(dtan(th0)*dcos(ph0)) < c thy0 = ratan(dtan(th0)*dsin(ph0)) < c --- > thindex = ng6indexx*xlam ! NG6 about 1.7mrad/ang > ph0 = 2.0d0*pi*vect(2) ! azimuthal angle (radians) > > if (thindex.lt.1E-7) thindex = 1E-7 > if (abs(1.d0-(1.d0-dcos(thindex))*vect(1)).gt.1.0d0) then > write(6,*) ' Arccos(theta_start)>1 ' > goto 321 > end if > > th0 = dacos(1.d0-(1.d0-dcos(thindex))*vect(1)) ! isotropic (0,thindex) bec > thx0 = ratan(dtan(th0)*dcos(ph0)) > thy0 = ratan(dtan(th0)*dsin(ph0)) > 815,816c832,833 < thx0 = 2*ng6indexx*xlam/pi*dasin(2.0d0*vect(1)-1.d0) < thy0 = 2*ng6indexy*xlam/pi*dasin(2.0d0*vect(2)-1.d0) --- > c thx0 = 2*ng6indexx*xlam/pi*dasin(2.0d0*vect(1)-1.d0) > c thy0 = 2*ng6indexy*xlam/pi*dasin(2.0d0*vect(2)-1.d0) 846c863,869 < call hf2(1110,real(x0),real(y0),1.) --- > ccccccc scale factor is fluxscale*totcnt from readimage.f *~1.028 since readimage active area > ccccccc from readimage (>10cunts) = 24.08cm2 turns out to be slightly smaller than > ccccccc active PSM area here 4*(2.75X2.25)....smaller effect than 24.75/24.08=1.028 --> 1.01 > cccccccc =1.2869E8*3.1356E5*1.01=4.076E12 > call hf2(2220,real(x0),real(y0),4.076E12/nevts) > if (dabs(y0).le.0.5) call hf1(441,real(x0),4.076E12/nevts*0.02) > if (dabs(x0).le.0.5) call hf1(442,real(y0),4.076E12/nevts*0.02) 895,898c918,922 < pl_del = 0.d0 < pltot = pl_del + pltot < thtot = gy*mag_B*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot < if (pltot.gt.500) write(6,*) 'Al 1 ',pltot --- > pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) > pltot = pl_del + pltot > thtot = gy*mag_B*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot > > if (pltot.gt.500) write(6,*) 'guide ',pltot 922c946 < if((x0*x1.lt.0d0).and.x1.gt.-1d5) --- > if((x0*x1.lt.0d0).and.x1.gt.-1d5) 929c953 < if((x0*x1.lt.0d0).and.x1.gt.-1d5) --- > if((x0*x1.lt.0d0).and.x1.gt.-1d5) 1006c1030 < pl_del = 0.d0 --- > pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) 1009c1033 < if (pltot.gt.500) write(6,*) 'Al 2 ',pltot --- > if (pltot.gt.500) write(6,*) 'Al ',pltot 1030c1054 < if (pltot.gt.500) write(6,*) 'Al input ',pltot --- > if (pltot.gt.500) write(6,*) 'guide ',pltot 1078,1080d1101 < INCOAVGflux = INCOAVGflux + 1/(ix2-ix1)/(iy2-iy1)/nevts < if ((2.0+x1)**2+y1**2.le.0.625**2) < & INCOPKflux = INCOPKflux + 1/pi/0.625**2/nevts 1083c1104 < call hf2(1111,real(x1),real(y1),1.) --- > 1132c1153 < if (pltot.gt.500) write(6,*) 'Al input ',pltot --- > if (pltot.gt.500) write(6,*) 'guide ',pltot 1140,1177c1161,1191 < ******************************************************************** < ******************************************************************** < * Neutrons passing through the gap (10.0 cm) between the input < * coil and the target. < * < c Update x,y,z,thx,thy,en values to new values < x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1 < c < call airgap(en0, x0, y0, z0, thx0, thy0, < & tairx1,tairx2,tairy1,tairy2,tairz2, en1, x1, y1, z1, < & thx1,thy1,airflg1) < < c call collimator(x0,x1,y0,y1) !Boroflex septum to prevent crossovers < < if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 < < call hf1(4,real(x1),1.) < call hf1(5017,real(180./3.14159*thx1),1.) < call hf1(5018,real(180./3.14159*thy1),1.) < if(thx1.ge.0..and.thy1.ge.0.) then < phi = ratan(dtan(thy1)/dtan(thx1)) < elseif(thx1.lt.0..and.thy1.ge.0.) then < phi = pi + ratan(dtan(thy1)/dtan(thx1)) < elseif(thx1.lt.0..and.thy1.lt.0.) then < phi = pi + ratan(dtan(thy1)/dtan(thx1)) < elseif(thx1.ge.0..and.thy1.lt.0.) then < phi = 2.*pi + ratan(dtan(thy1)/dtan(thx1)) < endif < c call hf1(5021,real(180./3.14159*phi),1.) < pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) < pltot = pl_del + pltot < thtot = gy*mag_B*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot < call hf1(6001,real(pltot),1.) < call hf1(6002,real(thtot),1.) < if (pltot.gt.500) write(6,*) 'gap 1 ',pltot < if (dabs(thtot).gt.1.0e9) then < write(6,*) ' thtot ', thtot, ', pl_del ',pl_del < write(6,*) ' en ',en1, ', pltot ',pltot --- > > ******************************************************************** > * Neutrons passing through the gap (1st 1/2 of 10.0 cm) between the input > * coil and the target. > * > c Update x,y,z,thx,thy,en values to new values > x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1 > c > call airgap(en0, x0, y0, z0, thx0, thy0, > & tairx1,tairx2,tairy1,tairy2,0.5*tairz2, en1, x1, y1, z1, > & thx1,thy1,airflg1) > > if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 > > > if(thx1.ge.0..and.thy1.ge.0.) then > phi = ratan(dtan(thy1)/dtan(thx1)) > elseif(thx1.lt.0..and.thy1.ge.0.) then > phi = pi + ratan(dtan(thy1)/dtan(thx1)) > elseif(thx1.lt.0..and.thy1.lt.0.) then > phi = pi + ratan(dtan(thy1)/dtan(thx1)) > elseif(thx1.ge.0..and.thy1.lt.0.) then > phi = 2.*pi + ratan(dtan(thy1)/dtan(thx1)) > endif > pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) > pltot = pl_del + pltot > thtot = gy*mag_B*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot > if (pltot.gt.500) write(6,*) 'gap 1 ',pltot > if (dabs(thtot).gt.1.0e9) then > write(6,*) ' thtot ', thtot, ', pl_del ',pl_del > write(6,*) ' en ',en1, ', pltot ',pltot 1179,1181c1193,1248 < end if < * < ******************************************************************** --- > end if > > ******************************************************************** > > ******************* INput Coil Exit Flux Measurement*************** > INCOAVGflux = INCOAVGflux + 1/(ix2-ix1)/(iy2-iy1)/nevts > if (x1**2+y1**2.le.0.625**2) > & INCOPKflux = INCOPKflux + 1/pi/0.625**2/nevts > ccccccc scale factor is fluxscale*totcnt from readimage.f *~1.028 since readimage active area > ccccccc from readimage (>10cunts) = 24.08cm2 turns out to be slightly smaller than > ccccccc active PSM area here 4*(2.75X2.25)....smaller effect than 24.75/24.08=1.028 --> 1.01 > cccccccc =1.2869E8*3.1356E5*1.01=4.076E12 > call hf2(2221,real(x1),real(y1),4.076E12/nevts) > if (dabs(y1).le.0.5) call hf1(443,real(x1),4.076E12/nevts*0.02) > if (dabs(x1).le.0.5) call hf1(444,real(y1),4.076E12/nevts*0.02) > ******************************************************************** > > ******************************************************************** > * Neutrons passing through the gap (2nd 1/2 of 10.0 cm) between the input > * coil and the target. > * > c Update x,y,z,thx,thy,en values to new values > x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1 > c > call airgap(en0, x0, y0, z0, thx0, thy0, > & tairx1,tairx2,tairy1,tairy2,0.5*tairz2, en1, x1, y1, z1, > & thx1,thy1,airflg1) > > if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 > > call hf1(4,real(x1),1.) > call hf1(5017,real(180./3.14159*thx1),1.) > call hf1(5018,real(180./3.14159*thy1),1.) > if(thx1.ge.0..and.thy1.ge.0.) then > phi = ratan(dtan(thy1)/dtan(thx1)) > elseif(thx1.lt.0..and.thy1.ge.0.) then > phi = pi + ratan(dtan(thy1)/dtan(thx1)) > elseif(thx1.lt.0..and.thy1.lt.0.) then > phi = pi + ratan(dtan(thy1)/dtan(thx1)) > elseif(thx1.ge.0..and.thy1.lt.0.) then > phi = 2.*pi + ratan(dtan(thy1)/dtan(thx1)) > endif > c call hf1(5021,real(180./3.14159*phi),1.) > pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) > pltot = pl_del + pltot > thtot = gy*mag_B*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot > call hf1(6001,real(pltot),1.) > call hf1(6002,real(thtot),1.) > if (pltot.gt.500) write(6,*) 'gap 1 ',pltot > if (dabs(thtot).gt.1.0e9) then > write(6,*) ' thtot ', thtot, ', pl_del ',pl_del > write(6,*) ' en ',en1, ', pltot ',pltot > write(6,*) ' z ', z1 > end if > > ******************************************************************** 1199c1266 < if (pltot.gt.500) write(6,*) 'Al 3 ',pltot --- > if (pltot.gt.500) write(6,*) 'guide ',pltot 1385c1452 < if (pltot.gt.500) write(6,*) 'Al 4 ',pltot --- > if (pltot.gt.500) write(6,*) 'guide ',pltot 1408c1475 < if (pltot.gt.500) write(6,*) 'Al Pi ',pltot --- > if (pltot.gt.500) write(6,*) 'guide ',pltot 1426a1494,1496 > pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) > pltot = pl_del + pltot > thtot = gy*mag_B*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot 1461c1531 < if (pltot.gt.500) write(6,*) 'Al Pi ',pltot --- > if (pltot.gt.500) write(6,*) 'guide ',pltot 1517c1587 < if (pltot.gt.500) write(6,*) 'Al 5 ',pltot --- > if (pltot.gt.500) write(6,*) 'guide ',pltot 1653,1654c1723,1728 < call hf1(2002,real(y1),wgt) < call hf2(2222,real(x1),real(y1),wgt) --- > call hf1(2002,real(y1),wgt) > ccccccc scale factor is fluxscale*totcnt from readimage.f *~1.028 since readimage active area > ccccccc from readimage (>10cunts) = 24.08cm2 turns out to be slightly smaller than > ccccccc active PSM area here 4*(2.75X2.25)....smaller effect than 24.75/24.08=1.028 --> 1.01 > cccccccc =1.2869E8*3.1356E5*1.01=4.076E12 > call hf2(2222,real(x1),real(y1),4.076E12*wgt/nevts) 1713c1787 < if (pltot.gt.500) write(6,*) 'Al 6 ',pltot --- > if (pltot.gt.500) write(6,*) 'guide ',pltot 1729,1732c1803 < & thx1,thy1,airflg1) < < c call collimator(x0,x1,y0,y1) !boroflex septum to prevent crossovers < --- > & thx1,thy1,airflg1) 1769c1840 < if (pltot.gt.500) write(6,*) 'Al output ',pltot --- > if (pltot.gt.500) write(6,*) 'guide ',pltot 1803,1806c1874,1876 < < if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 < call hf1(5,real(x1),wgt) < call hf2(3333,real(x1),real(y1),wgt) --- > > if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 > 1858c1928 < if (pltot.gt.500) write(6,*) 'Al output ',pltot --- > if (pltot.gt.500) write(6,*) 'Al ',pltot 1864a1935,1937 > call hf1(5,real(x1),wgt) > > 1875,1878c1948 < & thx1,thy1,airflg1) < < c call collimator(x0,x1,y0,y1) !boroflex septum to prevent crossovers < --- > & thx1,thy1,airflg1) 1885c1955 < if (x1**2+y1**2.le.0.625**2) --- > if ((1.5+x1)**2+y1**2.le.0.625**2) 1898c1968,1977 < --- > *************Output of OuputCoil Flux measurement****************** > ccccccc scale factor is fluxscale*totcnt from readimage.f *~1.028 since readimage active area > ccccccc from readimage (>10cunts) = 24.08cm2 turns out to be slightly smaller than > ccccccc active PSM area here 4*(2.75X2.25)....smaller effect than 24.75/24.08=1.028 --> 1.01 > cccccccc =1.2869E8*3.1356E5*1.01=4.076E12 > call hf2(3333,real(x1),real(y1),4.076E12*wgt/nevts) > if (dabs(y1).le.0.5) > & call hf1(445,real(x1),4.076E12*wgt/0.9976/nevts*0.02) > if (dabs(x1+0.9).le.0.5) > & call hf1(446,real(y1),4.076E12*wgt/0.9976/nevts*0.02) 1941c2020 < if (pltot.gt.500) write(6,*) 'Al 7 ',pltot --- > if (pltot.gt.500) write(6,*) 'guide ',pltot 1958,1961c2037 < < c call collimator(x0,x1,y0,y1) !boroflex septum to prevent crossovers < < if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 --- > if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 1967d2042 < if (pltot.gt.600) goto 100 2049c2124,2128 < call hf2(4444,real(x1),real(y1),wgt) --- > ccccccc scale factor is fluxscale*totcnt from readimage.f *~1.028 since readimage active area > ccccccc from readimage (>10cunts) = 24.08cm2 turns out to be slightly smaller than > ccccccc active PSM area here 4*(2.75X2.25)....smaller effect than 24.75/24.08=1.028 --> 1.01 > cccccccc =1.2869E8*3.1356E5*1.01=4.076E12 > call hf2(4444,real(x1),real(y1),4.076E12*wgt/nevts) 2579a2659 > 3393c3473 < write(6,*) ' num guide1 scatters ',gsc --- > write(6,*) ' num guide1 scatters ',gsc 3399c3479 < write(3,*) ' num guide1 scatters ',gsc --- > write(3,*) ' num guide1 scatters ',gsc