program main_WUp C Target in the Up position on the West side implicit none include 'const.inc' integer loop, dncross, nlamx, nlamy integer ngbx, ngbx0,ngbx1,ngbx2,ngbx3 integer ngby, ngby0,ngby1,ngby2,ngby3 integer nibx, nibx0,nibx1,nibx2,nibx3 integer niby, niby0,niby1,niby2,niby3 integer nobx, nobx0,nobx1,nobx2,nobx3 integer noby, noby0,noby1,noby2,noby3 integer nrfl1,nrfl2 integer scatflg1, scatflg2 double precision ep, wp, kp, k1, thp,php,qlast double precision dcrosw,dcrose,dscrosw,dscrose integer i, j, counter,lim,Event,cin,cot,trig !bec double precision Nenscat1,Ntar1,Nenscat2,Ntar2,Ntar12,NEdet,Nxo double precision encnt, Nsdet,Nthdet, Ntar1out, Ntar2out double precision Nalscat1, Nalscat2, Nthscat1, Nthscat2 double precision Nalout1, Nalout2, Nalin1, Nalin2, Nairin1 !rcm double precision Nairout1, Nairscat1 double precision NPSMout, Nincoilout, NOutcoilout, NTargetout ! bec double precision numtar1flg, numtar2flg, numendif1, numbeamsc double precision numendif2, numendifdet, numbounce double precision numbounce2, numbounce3,Ndetin double precision Nemp1in,Nemp2in,Nemp1out,Nemp2out, NtarWEin double precision INPKFlux,PSMAVGFlux,PSMPKFlux double precision INCOAVGFlux,INCOPKFlux,TARAVGFlux,TARPKFlux double precision OUTCOAVGFlux,OUTCOPKFlux,DETAVGFlux,DETPKFlux double precision xmax,ymax,capflpsm,capflinc,capflout,capfldet double precision fsw,fse,ufsw,ufse,rotnsw,rotnse,urotnsw,urotnse double precision rotsw,rotse,urotsw,urotse double precision fxw,fxe,ufxw,ufxe double precision rotxw,rotxe,urotxw,urotxe c real xlam, HRNDM1, vect(7), lm(3), const, xinit real yinit, plr, xtest, qvec(2) character*28 outfname double precision th0, ph0, thindex, enstart, entar1, entar2 double precision en0, x0, y0, z0, thx0, thy0, xstart double precision en1, x1, y1, z1, thx1, thy1 real xtmp, ytmp double precision time1 logical tar1flg, tar2flg,enflg1,enflg2,thflg1,thflg2,alflg1,alflg2 logical airflg1, histflg, histflg2, empflg1, empflg2, bounceflg logical bounceflg2, bounceflg3, beamscflg integer ranno integer seedc double precision Nempty1, Nempty2 C parameter(ranno=1) double precision r2d, lambda, sig, intlen, thscat double precision phi, ztar1,scwgt double precision pltot,thtot,pl_del,rotang double precision lsrot,lsrotsq,lsrotavg,lsrottsc,lsrottscsq double precision lsrotgsc,lsrotgscsq,lsde,lsdesq double precision lsrotxsc,lsrotxscsq double precision rsrot,rsrotsq,rsrotavg,rsrottsc,rsrottscsq double precision rsrotnsc,rsrotnscsq,lsrotnsc,lsrotnscsq double precision rsrotgsc,rsrotgscsq,rsde,rsdesq double precision rsrotxsc,rsrotxscsq double precision thxinput, thyinput double precision nr,nrts,nrgs,nl,nlts,nlgs,nlxs,nrxs,nder,ndel,nb double precision nrsq,nrtssq,nrgssq,nlsq,nltssq,nlgssq double precision nrxssq,nlxssq double precision ndersq,ndelsq,nrtsraw,nltsraw,nrgsraw,nlgsraw double precision nrxsraw,nlxsraw double precision nrnsc,nlnsc,numE2W, numW2E,nrnscsq,nlnscsq double precision lall,rall,lallsq,rallsq,lsc,rsc,lscsq,rscsq double precision nrPA,nrPAsq,rPAavg,nlPA,nlPAsq,lPAavg !rcm PAsin double precision uncrsrotavg,uncrPAavg,unclsrotavg,unclPAavg double precision nrnscPA,nrnscPAsq,rsnscPAavg,rsrotnscavg double precision nlnscPA,nlnscPasq,lsnscPAavg,lsrotnscavg double precision uncrsnsc,uncrsnscPA,unclsnsc,unclsnscPA double precision nrtsPA,nrtsPAsq,rtsPAavg,rsrottscavg double precision nltsPA,nltsPAsq,ltsPAavg,lsrottscavg double precision uncrstsc,uncrstscPA,unclstsc,unclstscPA double precision nrxsPA,nrxsPAsq,rxsPAavg,rsrotxscavg double precision nlxsPA,nlxsPAsq,lxsPAavg,lsrotxscavg double precision uncrsxsc,uncrsxscPA,unclsxsc,unclsxscPA double precision nrgsPA,nrgsPAsq,rgsPAavg,rsrotgscavg double precision nlgsPA,nlgsPAsq,lgsPAavg,lsrotgscavg double precision uncrsgsc,uncrsgscPA,unclsgsc,unclsgscPA double precision nrallPA,nrallPAsq,nlallPA,nlallPAsq,rallavg double precision rallPAavg,lallavg,lallPAavg,nrscPA,rscPAavg double precision rscavg,nlscPA,lscPAavg,lscavg,nlscPAsq,nrscPAsq double precision uncrallavg,uncrallPAavg,unclallavg,unclallPAavg double precision argrsrot,argrsrotPA,argrsnsc,argrsnscPA double precision uncrscavg,uncrscPAavg,unclscavg,unclscPAavg double precision arglsrot,arglsrotPA,arglsnsc,arglsnscPA double precision argltsc,argltscPA,arglgsc,arglgscPA double precision argrtsc,argrtscPA,argrgsc,argrgscPA double precision arglxsc,arglxscPA double precision argrxsc,argrxscPA double precision argrall,argrallPA,arglall,arglallPA double precision zdet, rsqrt,ratan,Nguideout double precision NguideoutE,NguideoutW,NincoiloutE,NincoiloutW double precision NoutcoiloutW,NoutcoiloutE double precision zlm integer beamflg,nbx,nby,scflag,gsc,isc,osc logical nanflg real beamlength c---- double precision lwde,rwde,lwdesq,rwdesq,lpl,lplsq,rpl,rplsq double precision lplns,rplns,lplnssq,rplnssq integer nlall,nrall,nlsc,nrsc,ndelw,nderw c---- double precision IntBz,temp,norm !for gradient runs -KG double precision PA, xtar1,xtar2 !Polarization product -KG CHARACTER day*2,month*2,hour*2,minute*2, filename*20 c----------------------------------- integer ispec double precision spec(1000,2), alambda c----------------------------------- integer qmin, qmax, hmax integer hws1, hws2 double precision k0, qdummy, q, q1, wp1, ep1, kp1, th_test1 double precision q2, wp2, ep2, kp2, th_test2 double precision qtest, wtest real wgt, wgt2 c----------------------------------- c c ------------------------- c fill the ntuple 'spin.rz' //// ntuple variables .... integer max_paw, iquest integer istat, icycle, datetime(8) parameter(max_paw = 2000000) ! PAW common block size c parameter(max_paw = 1000000) ! PAW common block size real quest, p character date*12, time*12, zone*12 c common/quest/iquest(100) common/pawc/p(max_paw) common/bfield/temp(1101),norm(1101) if (bzgradnum.eq.1) then open(unit=37,file=bzfile1,status='old') else if (bzgradnum.eq.2) then open(unit=37,file=bzfile2,status='old') else if (bzgradnum.eq.3) then open(unit=37,file=bzfile3,status='old') else if (bzgradnum.eq.4) then open(unit=37,file=bzfile4,status='old') end if if (bzgradnum.ne.0) then read(37,*) do nb = 1, 1101 read(37,*) temp(nb) c write(6,*) temp(nb) end do read(37,*) do nb = 1, 1101 read(37,*) norm(nb) c write(6,*) norm(nb) end do end if close(unit=37) c time1=time() c ranno=mod(time1,2000000) c ranno=mod(time1,1089926518) CALL Date_and_time(date,time,zone,datetime) !bec day= CHAR(INT(datetime(3)/10)+48)//CHAR(mod(datetime(3),10)+48) month= CHAR(int(datetime(2)/10)+48)//CHAR(MOD(datetime(2),10)+48) hour= CHAR(INT(datetime(5)/10)+48)//CHAR(mod(datetime(5),10)+48) minute= CHAR(INT(datetime(6)/10)+48)//CHAR(mod(datetime(6),10)+48) filename ='xportWUp-'//month//day//'-'//hour//minute//'.o' open(unit=3,FILE=filename,Status='new') write(*,*) ' Nevents: ', nevts ! max_paw !bec write(*,*) ' Start time: ', date,time !bec write(3,*) ' Nevents: ', nevts ! max_paw !bec write(3,*) ' Start time: ', date,time !bec write(*,*) ' Upstream target on West (Right), x>0.' write(3,*) ' Upstream target on West (Right), x>0.' c print *,"asking for seed" c read(*,*) ranno; c ranno=seedc() ranno=datetime(6)*datetime(7)*datetime(8)+datetime(7)+datetime(8) !bec write(*,*) 'ranNo is ', ranno write(3,*) 'ranNo is ', ranno c stop c stop/N/B/msarsour/spin/new_ih/ call hlimit(max_paw) c iquest(10) = 6500000 c c--------------------------- beamlength = talz1+talz2+2.0*talz3+2.0*talz4+2.0*talz5+ & 4.0*talz6+talz7+ & tairz1+tairz2+tairz3+tairz4+tairz5+tairz6+ & gl+2.0*tarz+iz0+oz+sz write(6,*) ' Calculated beamlength = ',beamlength c----------------------------------------------------------------------- c----------------------------------------------------------------------- c call hrget(99,'hbooks/lam.hbook','') !New initial spectrum after PSM -KG c call hrget(102,'hbooks/xyspot.hbook','') call hrget(151,'hbooks/gradx.hbook','') call hrget(152,'hbooks/gradx.hbook','') call hrget(161,'hbooks/ng6-parab-xy.hbook','') call hrget(162,'hbooks/ng6-parab-xy.hbook','') c call hrget(161,'hbooks/ngc-parab8x8-xy.hbook','') c call hrget(162,'hbooks/ngc-parab8x8-xy.hbook','') c call hrget(161,'hbooks/ngc-parab11x11-xy.hbook','') c call hrget(162,'hbooks/ngc-parab11x11-xy.hbook','') call hrget(699,'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 cc----------------------------------------------------------------------- c read(*,*) ranno call rluxgo(3,ranno,0,0) ! RANLUX luxury level 2 c----------------------------------------------------------------------- c Print *, "setting up Hbooks" call hbook1(101,'lamda_start',100,0,16.0,0.) !-KG call hbook1(305,'West q sim(1/Ang)',100,0,0.1,0.) ! bec test q choices call hbook1(306,'West w sim(1/psec)',100,-0.5,0.5,0.) ! bec test q choices call hbook1(307,'West Eo sim(meV)',3000,0.,40.0,0.) ! bec test q choices call hbook1(308,'West Ep sim(meV)',3000,0.,40.0,0.) ! bec test q choices call hbook1(309,'West ko sim(1/Ang)',100,0,4.0,0.) ! bec test q choices call hbook1(310,'West kp sim(1/Ang)',100,0,4.0,0.) ! bec test q choices call hbook1(205,'East q sim(1/Ang)',100,0,0.1,0.) ! bec test q choices call hbook1(206,'East w sim(1/psec)',100,-0.5,0.5,0.) ! bec test q choices call hbook1(207,'East Eo sim(meV)',3000,0.,40.0,0.) ! bec test q choices call hbook1(208,'East Ep sim(meV)',3000,0.,40.0,0.) ! bec test q choices call hbook1(209,'East ko sim(1/Ang)',100,0,4.0,0.) ! bec test q choices call hbook1(210,'East kp sim(1/Ang)',100,0,4.0,0.) ! bec test q choices call hbook1(1001,'x (cm) Uptar ',100,-5.5,5.5,0.) call hbook1(1,'x (cm) start',100,-5.5,5.5,0.) call hbook1(10,'y (cm) start',100,-5.5,5.5,0.) call hbook1(2,'x (cm) guide out',100,-5.5,5.5,0.) call hbook1(3,'x (cm) incoil out',100,-5.5,5.5,0.) call hbook1(4,'x (cm) target in',100,-5.5,5.5,0.) call hbook1(5,'x (cm) outcoil out',100,-5.5,5.5,0.) call hbook1(6,'x (cm) ASM out',100,-5.5,5.5,0.) call hbook1(20,'E (meV) start Both',2000,0.,40.,0.) call hbook1(31,' thx start (rad)',200,-0.1,0.1,0.) call hbook1(32,' thy start (rad)',200,-0.1,0.1,0.) 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.) call hbook1(1002,'y (cm) UpWtar',100,-6.,6.,0.) c call hbook1(1003,'thx (mrad) Up',160,-1.6,1.6,0.) c call hbook1(1004,'thy (mrad) Up',160,-1.6,1.6,0.) call hbook1(1005,'E (meV)/ non-scat tar1',2000,0.,40.,0.) c call hbook1(1006,'Path len (cm)/ non-scat UpB',100,35.0,100.0,0.) c call hbook1(1007,'th_rot (mrad)/ non-scat UpB',100,0.,10000.,0.) c call hbook1(1008,'Path len (cm)/ non-scat UpR',100,35.0,100.,0.) c call hbook1(1009,'th_rot (mrad)/ non-scat UpR',100,0.,10000.,0.) c call hbook1(1010,'Path len (cm)/ non-scat UL',100,35.0,100.0,0.) c call hbook1(1011,'th_rot (mrad)/ non-scat UL',100,0.,10000.,0.) call hbook1(1012,'E (meV)/ scat UpWE',2000,0.,40.,0.) c call hbook1(1013,'Path len (cm)/ scat UpB',100,35.0,100.0,0.) c call hbook1(1014,'th_rot (mrad)/ scat UpB',100,0.,10000.,0.) c call hbook1(1015,'Path len (cm)/ scat UpR',100,35.0,100.0,0.) c call hbook1(1016,'th_rot (mrad)/ scat UpR',100,0.,10000.,0.) c call hbook1(1017,'Path len (cm)/ scat UpL',100,35.0,100.0,0.) c call hbook1(1018,'th_rot (mrad)/ scat UpL',100,0.,10000.,0.) * call hbook1(2001,'x (cm) Dn',100,-5.5,5.5,0.) call hbook1(2002,'y (cm) Dn',100,-6.,6.,0.) c call hbook1(2003,'thx (mrad) Dn',160,-1.6,1.6,0.) c call hbook1(2004,'thy (mrad) Dn',160,-1.6,1.6,0.) call hbook1(2005,'E (meV)/ non-scat DnWE',2000,0.,40.,0.) c call hbook1(2006,'Path len (cm)/ non-scat DnB',100,41.5,41.7,0.) c call hbook1(2007,'th_rot (mrad)/ non-scat DnB',100,0.,1000.,0.) c call hbook1(2008,'Path len (cm)/ non-scat DnR',100,41.5,41.7,0.) c call hbook1(2009,'th_rot (mrad)/ non-scat DnR',100,0.,1000.,0.) c call hbook1(2010,'Path len (cm)/ non-scat DnL',100,41.5,41.7,0.) c call hbook1(2011,'th_rot (mrad)/ non-scat DnL',100,0.,1000.,0.) call hbook1(2012,'E (meV)/ scat DnWE',2000,0.,40.,0.) c call hbook1(2013,'Path len (cm)/ scat DnB',100,41.,45.,0.) c call hbook1(2014,'th_rot (mrad)/ scat DnB',100,0.,1000.,0.) c call hbook1(2015,'Path len (cm)/ scat DnR',100,41.,45.,0.) c call hbook1(2016,'th_rot (mrad)/ scat DnR',100,0.,1000.,0.) c call hbook1(2017,'Path len (cm)/ scat DnL',100,41.,45.,0.) c call hbook1(2018,'th_rot (mrad)/ scat DnL',100,0.,1000.,0.) * call hbook1(3001,'x (cm) Det',100,-3.5,3.5,0.) call hbook1(3002,'y (cm) Det',100,-6.,6.,0.) call hbook1(3003,'thx (rad) Det',200,-0.1,0.1,0.) call hbook1(3004,'thy (rad) Det',200,-0.1,0.1,0.) call hbook1(3005,'E (meV) DetWE',2000,0.,40.,0.) call hbook1(4004,'Lambda (Angs) Initial',1000,0.,16.,0.) !KG call hbook1(4005,'Lambda (Angs) DetWE',1000,0.,16.,0.) c call hbook1(3006,'Path length (cm) DeWE',200,480.58,480.70,0.) !KG c call hbook1(3007,'th_rotation (mrad) DetWE',200,-400.,0.0,0.) !KG c call hbook1(3008,'Path length (cm) DetEast',200,114.58,114.7,0.) c call hbook1(3008,'Path length (cm) DetEast',200,480.50,480.81,0.) c call hbook1(3009,'th_rotation (mrad) DetEast',200,-5.,5.0,0.) c call hbook1(3010,'Path length (cm) DetWest',200,114.58,114.7,0.) c call hbook1(3010,'Path length (cm) DetWest',200,480.50,480.81,0.) c call hbook1(3011,'th_rotation (mrad) DetWest',200,-5.,5.0,0.) call hbook1(3012,'E (meV) DetEast',2000,0.,40.,0.) call hbook1(3013,'E (meV) DetWest',2000,0.,40.,0.) c call hbook1(3014,'E (meV) scat DetWE',2000,0.,40.,0.) c call hbook1(3020,'E_loss (meV) scat DetWE',2000,-0.4,0.4,0.) c call hidopt(3020,'logy') call hbook1(3021,'E_loss (meV) scat DetE',2000,-0.5,0.5,0.) call hbook1(3022,'E_loss (meV) scat DetW',2000,-0.5,0.5,0.) c c adding more spectra below for diagnostics bec call hbook2(8010,'sig tot (b/atom) vs. q (1/Ang)',50,0.0,0.2, & 50,0.0,1.0,0.) call hbook1(4010,'West sig tot (b/atom)',200,0.,1.,0.) c call hbook1(4011,'West interaction len (cm)',200,0.,2000.,0.) call hbook1(4020,'East sig tot (b/atom)',200,0.,1.,0.) c call hbook1(4021,'East interaction len (cm)',200,0.,2000.,0.) call hbook1(4012,'theta scattered (deg) U',200,0.,180.,0.) c call hbook1(4015,'theta-p (deg) U',200,0.,180.,0.) c call hbook1(4016,'phi-p (deg) U',200,0.,360.,0.) call hbook1(4017,'th-0 (deg)',200,0.,10.,0.) c call hbook1(4018,'phi-0 (deg) ',200,0.,360.,0.) call hbook1(5012,'theta scattered (deg) D',200,0.,180.,0.) c call hbook1(5015,'theta-p (deg) D',200,0.,180.,0.) c call hbook1(5016,'phi-p (deg) D',200,0.,360.,0.) call hbook1(5017,'theta-x (deg) U',200,-2.,2.,0.) call hbook1(5018,'theta-y (deg) U',200,-2.,2.,0.) call hbook1(5019,'theta-x (deg) D',200,-2.,2.,0.) call hbook1(5020,'theta-y (deg) D',200,-2.,2.,0.) c call hbook1(5021,'phi (deg) af-guide',200,0.,360.,0.) c call hbook1(5022,'phi (deg) af-target 1',200,0.,360.,0.) c call hbook1(5023,'phi (deg) af-target 2',200,0.,360.,0.) call hbook1(6001,'Path length (cm) enter Up WE',200,212.0, & 215.0,0.) c call hbook1(6003,'Path length (cm) Det scat WE',300,480.50, c & 480.90,0.) c call hbook1(6004,'th_rotation (mrad) Det scat WE',400,-800., c & 0.0,0.) call hbook1(6005,'Path length (cm) Det scat West',300, & beamlength-0.5,beamlength+2.0,0.) call hbook1(6007,'Path length (cm) Det scat East',300, & beamlength-0.5,beamlength+2.0,0.) c call hbook1(6009,'Path length (cm) Det no scat WE',300,480.50, c & 480.90,0.) c call hbook1(6010,'th_rotation (mrad) Det no scat WE',400,-800., c & 0.0,0.) call hbook1(6011,'Path length (cm) Det no scat East',300, & beamlength-0.5,beamlength+2.0,0.) call hbook1(6013,'Path length (cm) Det no scat West',300, & beamlength-0.5,beamlength+2.0,0.) call hbook1(6015,'Crossover Path length (cm) Det scat West', & 300,beamlength-0.5,beamlength+2.0,0.) call hbook1(6017,'Crossover Path length (cm) Det scat East', & 300,beamlength-0.5,beamlength+2.0,0.) call hbook1(6035,'Crossover Path length (cm) DetEast', & 300,beamlength-0.5,beamlength+2.0,0.) call hbook1(6037,'Crossover Path length (cm) DetWest', & 300,beamlength-0.5,beamlength+2.0,0.) c c Bz gradient and Pi-coil adjustments c if (bzgradnum.ne.0) then c if (picoilflg) then ! Bz gradient, Picoil ON c call hbook1(3009,'th_rotation (mrad) DetEast',300,-800.,10.0,0.) c call hbook1(3011,'th_rotation (mrad) DetWest',300,-800.,10.0,0.) call hbook1(6002,'th_rotation (mrad) enter Up WE',400,-100., & 100.,0.) call hbook1(6006,'th_rotation (mrad) Det scat West',400,-800., & 0.0,0.) call hbook1(6008,'th_rotation (mrad) Det scat East',400,-800., & 0.0,0.) call hbook1(6012,'th_rotation (mrad) Det no scat East',400,-800., & 0.0,0.) call hbook1(6014,'th_rotation (mrad) Det no scat West', & 400,-800.,0.,0.) call hbook1(6016,'Crossover th_rotation (mrad) Det scat West', & 400,-800.,0.0,0.) call hbook1(6018,'Crossover th_rotation (mrad) Det scat East', & 400,-800.,0.0,0.) call hbook1(6036,'Crossover th_rotation (mrad) DetEast', & 400,-800.,0.0,0.) call hbook1(6038,'Crossover th_rotation (mrad) DetWest', & 400,-800.,0.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.) else ! Bz gradient, Picoil OFF c call hbook1(3009,'th_rotation (mrad) DetEast',300,-10.,80.0,0.) c call hbook1(3011,'th_rotation (mrad) DetWest',300,-10.,80.0,0.) call hbook1(6002,'th_rotation (mrad) enter Up WE',400,0., & 60.,0.) call hbook1(6006,'th_rotation (mrad) Det scat West',400,0., & 60.0,0.) call hbook1(6008,'th_rotation (mrad) Det scat East',400,0., & 60.0,0.) call hbook1(6012,'th_rotation (mrad) Det no scat East',400,0., & 60.0,0.) call hbook1(6014,'th_rotation (mrad) Det no scat West', & 400,0.,60.,0.) call hbook1(6016,'Crossover th_rotation (mrad) Det scat West', & 400,0.,60.0,0.) call hbook1(6018,'Crossover th_rotation (mrad) Det scat East', & 400,0.,60.0,0.) call hbook1(6036,'Crossover th_rotation (mrad) DetEast', & 400,0.,60.0,0.) call hbook1(6038,'Crossover th_rotation (mrad) DetWest', & 400,0.,60.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,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.) end if else ! No Bz gradient If (picoilflg) then ! NO Bz gradient, Picoil ON c call hbook1(3009,'th_rotation (mrad) DetEast',300,-10.,10.0,0.) call hbook1(3011,'th_rotation (mrad) DetWest',300,-10.,10.0,0.) call hbook1(6002,'th_rotation (mrad) enter Up WE',300,-10., & 20.,0.) call hbook1(6006,'th_rotation (mrad) Det scat West',300,-10., & 10.0,0.) call hbook1(6008,'th_rotation (mrad) Det scat East',300,-10., & 10.0,0.) call hbook1(6012,'th_rotation (mrad) Det no scat East',300,-10., & 10.0,0.) call hbook1(6014,'th_rotation (mrad) Det no scat West', & 300,-10.,10.,0.) call hbook1(6016,'Crossover th_rotation (mrad) Det scat West', & 300,-10.,10.0,0.) call hbook1(6018,'Crossover th_rotation (mrad) Det scat East', & 300,-10.,10.0,0.) call hbook1(6036,'Crossover th_rotation (mrad) DetEast', & 300,-10.,10.0,0.) call hbook1(6038,'Crossover th_rotation (mrad) DetWest', & 300,-10.,10.0,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.) else ! NO Bz gradient, Picoil OFF call hbook1(3009,'th_rotation (mrad) DetEast',300,-5.,50.0,0.) call hbook1(3011,'th_rotation (mrad) DetWest',300,-5.,50.0,0.) call hbook1(6002,'th_rotation (mrad) enter Up WE',300,-5., & 50.,0.) call hbook1(6006,'th_rotation (mrad) Det scat West',300,-5., & 50.0,0.) call hbook1(6008,'th_rotation (mrad) Det scat East',300,-5., & 50.0,0.) call hbook1(6012,'th_rotation (mrad) Det no scat East',300,-5., & 50.0,0.) call hbook1(6014,'th_rotation (mrad) Det no scat West', & 300,-5.,50.,0.) call hbook1(6016,'Crossover th_rotation (mrad) Det scat West', & 300,-5.,50.0,0.) call hbook1(6018,'Crossover th_rotation (mrad) Det scat East', & 300,-5.,50.0,0.) call hbook1(6036,'Crossover th_rotation (mrad) DetEast', & 300,-5.,50.0,0.) call hbook1(6038,'Crossover th_rotation (mrad) DetWest', & 300,-5.,50.0,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.) end if end if call hbook1(7001,'dE / 1st target',400,-1.,1.,0.) call hbook1(7002,'dE / 2nd target',400,-1.,1.,0.) call hbook1(7005,' Wgt Factor All',200,0.,2.,0.) call hbook1(7006,' Wgt Factor Any Scat',200,0.,1.1,0.) call hbook1(7007,' Wgt Factor Tar Scat',200,0.,1.1,0.) c call hbook2(8011,'Thx vs. E (meV) scat DetB',100,0.,40., c & 200,-0.1,0.1,0.) c call hbook2(8012,'Thx vs. E-loss (meV) scat DetB',100, c & -0.4,0.4,200,-0.1,0.1,0.) call hbook2(9999,'PA vs.WaveLength (Ang)',160,0.,16.,102, & 0.,1.02,0.) !-KG c call hbook2(9998,'omega (1/psec) vs. q (1/Ang)',100,0.,0.2,100, c & -2.0,2.0,0.) !-KG 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 Nenscat1=0 !bec Ntar1=0 !bec Nthscat1=0 Nenscat2=0 !bec Ntar2=0 !bec Nthscat2=0 Nsdet=0 !bec NEdet=0 !bec nrtsraw=0;nltsraw=0;nrgsraw=0;nlgsraw=0 NTargetout=0; Ntar1out=0; Ntar2out=0 !bec Nthdet=0 NPSMout=0 ; Nincoilout=0 ; NOutcoilout = 0 !bec lsrot=0;lsrotsq=0;lsrottsc=0;lsrottscsq=0;lsrotgsc=0;lsrotgscsq=0 lsde=0;lsdesq=0;rsrot=0;rsrotsq=0;rsrottsc=0;rsrottscsq=0; rsrotavg=0;lsrotavg=0 nrnsc=0;nlnsc=0;rsrotnsc=0;lsrotnsc=0;rsrotnscsq=0;lsrotnscsq=0 rsrotgsc=0 rsrotgscsq=0;rsde=0;rsdesq=0;nr=0;nrts=0;nrgs=0;nl=0;nlts=0;nlgs=0 nrPA=0;nrPAsq=0;rPAavg=0;nlPA=0;nlPAsq=0;lPAavg=0 !rcm PAsin uncrsrotavg=0;uncrPAavg=0;unclsrotavg=0;unclPAavg=0 nrnscPA=0;nrnscPAsq=0;rsnscPAavg=0;rsrotnscavg=0 nlnscPA=0;nrnscPAsq=0;lsnscPAavg=0;lsrotnscavg=0 uncrsnsc=0;uncrsnscPA=0;unclsnsc=0;unclsnscPA=0 nrtsPA=0;nrtsPAsq=0;rtsPAavg=0;rsrottscavg=0;uncrstsc=0 uncrstscPA=0;nltsPA=0;nltsPAsq=0;ltsPAavg=0;lsrottscavg=0 uncrstsc=0;uncrstscPA=0;unclstsc=0;unclstscPA=0 nrgsPA=0;nrgsPAsq=0;rgsPAavg=0;rsrotgscavg=0;uncrsgsc=0; uncrsgscPA=0;nlgsPA=0;nlgsPAsq=0;lgsPAavg=0;lsrotgscavg=0 uncrsxscPA=0; unclsxscPA=0 unclsgsc=0;unclsgscPA=0;nrallPA=0;nrallPAsq=0;nlallPA=0;rallavg=0 uncrscavg=0;uncrscPAavg=0;unclscavg=0;unclscPAavg=0 rallPAavg=0;lallavg=0;lallPAavg=0;nrscPA=0;rscPAavg=0;rscavg=0 nlscPA=0;lscPAavg=0;lscavg=0;nlallPAsq=0;uncrallavg=0; uncrallPAavg=0;unclallavg=0;unclallPAavg=0;nlscPAsq=0;nrscPAsq=0 nder=0;ndel=0 Nalscat1=0;Nalscat2=0;Nalout1=0;Nalout2=0; Nalin1=0; Nalin2=0 !rcm tar2flg = .false. ! bec Nairin1 = 0; Nairout1 = 0; Nairscat1 = 0 histflg = .false. ; histflg2 = .false. numtar1flg = 0; numendif1 = 0 numtar2flg = 0; numendif2 = 0 ; numbeamsc = 0.0 Nemp1in=0; Nemp2in=0; Nemp1out=0; Nemp2out=0; NtarWEin=0 Nxo=0 !-KG cross overs numE2W=0 ; numW2E=0 ;numbounce=0.0 lwde=0.0;rwde=0.0;lwdesq=0.0;rwdesq=0.0 nlall=0;nrall=0;nlsc=0;nrsc=0;ndelw=0;nderw=0 ngbx=0;ngbx0=0;ngbx1=0;ngbx2=0;ngbx3=0 ngby=0;ngby0=0;ngby1=0;ngby2=0;ngby3=0 xmax = 0.d0; ymax=0.d0 INPKFlux=0.d0;PSMAVGFlux=0.d0;PSMPKFlux=0.d0;INCOAVGFlux=0.d0 INCOPKFlux=0.d0;TARAVGFlux=0.d0;TARPKFlux=0.d0;OUTCOAVGFlux=0.d0 OUTCOPKFlux=0.d0;DETAVGFlux=0.d0;DETPKFlux=0.d0 nanflg=.false. ; scwgt = 1.0d0 c do 101 i = 1, nevts Event = 1 pltot = 0.0 thtot = 0.0 enflg1=.false. thflg1=.false. enflg2=.false. thflg2=.false. wgt = 1.0 ! assume non-scattered until see scattering in a target wgt2 = 1.0 ! assume non-scattered until see scattering in a target tar1flg = .false. ; tar2flg = .false. ; bounceflg = .false. bounceflg2 = .false. ; bounceflg3 = .false. scatflg1 = 0 ; scatflg2 = 0 ; beamscflg = .false. c c Generate initial neutron wavelength c c if (lamflg) then if (i.eq.1) then write(*,*) ' Lambda from Distribution ' write(3,*) ' Lambda from Distribution ' endif xlam = HRNDM1(99) !New initial spectrum after PSM -KG call hf1(101,real(xlam),1.) !-KG if(xlam.lt.1.) goto 100 if(xlam.gt.20) goto 100 else if (i.eq.1) then write(*,*) ' All Lambda = 5Ang ' write(3,*) ' All Lambda = 5Ang ' endif xlam = 5.0d0 end if c c call hf1(4004,real(xlam),1.) !Initial lambda histogram -KG c cbec c c The energy calculated from the neutron's wavelength c en0 = (1/(2*m_n))*(2*pi*h_par/xlam)**2 enstart = en0 if (en0.lt.0.d0) then write(6,*) ' en<0, en=1e-6' goto 100 end if call hf1(20,real(en0),1.) if (en0.le.1.0) encnt =encnt+1 c c Generates the initial trajectory of the neutron. c 321 call ranlux(vect,7) c if (vect(5).eq.0.50) write(*,*) 'vect(5)=0.50' !for wq.f -KG c vect(6) determines initial E or W, vect(7) gives the exact postition -KG c c thindex = ng6index*xlam ! NG6 about 1.7mrad/ang 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 c th0 = dacos(1.d0-(1.d0-dcos(thindex))*vect(1)) ! isotropic (0,thindex) bec cc write(99,*) vect(1), th0 c c if (th0.lt.1.E-8) th0=1.E-8 c c ph0 = 2.0d0*pi*vect(2) ! azimuthal angle (radians) c call hf1(4017,real(th0*180./3.14159),1.) cc call hf1(4018,real(ph0*180./3.14159),1.) c if (i.eq.1) then if (hegastar) write(6,*) ' !! He gas target run !! ' if (hegastar) write(3,*) ' !! He gas target run !! ' if (airtar) write(6,*) ' !! Air target run !! ' if (airtar) write(3,*) ' !! Air target run !! ' if (bzgradnum.eq.0) then write(6,*) ' ! NO Bz gradient in target region !' write(3,*) ' ! NO Bz gradient in target region !' else if (bzgradnum.eq.1) then write(6,*) ' !!!!Running with Bz gradient: ',bzfile1,' !!!!' write(3,*) ' !!!!Running with Bz gradient: ',bzfile1,' !!!!' else if (bzgradnum.eq.2) then write(6,*) ' !!!!Running with Bz gradient: ',bzfile2,' !!!!' write(3,*) ' !!!!Running with Bz gradient: ',bzfile2,' !!!!' else if (bzgradnum.eq.3) then write(6,*) ' !!!!Running with Bz gradient: ',bzfile3,' !!!!' write(3,*) ' !!!!Running with Bz gradient: ',bzfile3,' !!!!' else if (bzgradnum.eq.4) then write(6,*) ' !!!!Running with Bz gradient: ',bzfile4,' !!!!' write(3,*) ' !!!!Running with Bz gradient: ',bzfile4,' !!!!' end if if (picoilflg) then write(6,*) ' ! pi-coil ON ! ' write(3,*) ' ! pi-coil ON ! ' else write(6,*) ' ! pi-coil OFF ! ' write(3,*) ' ! pi-coil OFF ! ' end if if (sptmflg) then write(6,*) ' ! Septum Upstream ! ' write(3,*) ' ! Septum Upstream ! ' else write(6,*) ' ! No Septum Upstream ! ' write(3,*) ' ! No Septum Upstream ! ' end if write(6,*) write(3,*) end if ccccccccccccccccccccccc No x-y gradient if (i.eq.1) then write(*,*) ' ! Running with Flat xy-distribution !' write(3,*) ' ! Running with Flat xy-distribution !' end if xtest = xin*(2.*vect(3)-1.) !-KG x0 = xtest xstart = x0 !-KG y0 = yin*(2.*vect(4)-1.) ccccccccc gradient cc if (i.eq.1) then cc write(*,*) ' !!! Running with Position (x) gradient !' cc write(3,*) ' !!! Running with Position (x) gradient !' cc end if 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 c if (i.eq.1) then c write(*,*) ' ! xy-distribution from psm.hbook !' c write(3,*) ' ! xy-distribution from psm.hbook !' c end if c call HRNDM2(200,xtmp,ytmp) c x0= dble(xtmp) c y0 = dble(ytmp) c xtest = x0 c xstart = x0 c c c if ((abs(x0).gt.xin).or.(abs(y0).gt.yin)) goto 100 c c if ((x0**2+y0**2).le.0.625**2) c & INPKflux = INPKflux + 1/pi/0.625**2/nevts c c c cccccccccccc ccccccccc Inverted Parabola x - y distributions c 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 c c if (x0**2+y0**2.le.0.625**2) c & INPKflux = INPKflux + 1/pi/0.625**2/nevts c c ccccccccc z0= 0.0d0 cc cc Calculate the incident angles in the xz-plane and yz-plane cc note: these are not "direction cosines"= angle between momentum and x (or y) cc c thx0 = ratan(dtan(th0)*dcos(ph0)) c thy0 = ratan(dtan(th0)*dsin(ph0)) CC************************************************************************ c************for nSpinSim3.1.3 choosing thx0 and thy0 evening withing their c respective critical angles.. c isotropic c c if (i.eq.1) then c write(6,*) ' ! Isotropic angles!',ng6indexx,ng6indexy c write(3,*) ' ! Isotropic angles!',ng6indexx,ng6indexy c end if c c 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 c flat c c if (i.eq.1) then c write(6,*) ' ! Flat angular dist.!',ng6indexx,ng6indexy c write(3,*) ' ! Flat angular dist.!',ng6indexx,ng6indexy c end if c thx0 = ng6indexx*xlam*(vect(1)-0.5d0)*2.0d0 c thy0 = ng6indexy*xlam*(vect(2)-0.5d0)*2.0d0 c c cosine dist if (i.eq.1) then write(6,*); write(3,*) write(6,*) ' ! Cosine Angular Dist.!',ng6indexx,ng6indexy write(3,*) ' ! Cosine Angular Dist.!',ng6indexx,ng6indexy write(6,*); write(3,*) end if !!!!set maximum divergence in const.inc ng6indexx,y 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 C********for nSPinSim3.1.4 choosing thx0 and thy0 from hbooks, interpolated C******distributions from J. Cook simulations of NG6 and NGC c c if (i.eq.1) then c write(6,*) ' ! Angular Dist from J. Cook !' c write(3,*) ' ! Angular Dist from J. Cook !' c end if c if(xlam.lt.2.or.xlam.gt.12.) goto 101 c nlamx = int(xlam*10+1)+980 c nlamy = int(xlam*10+1)+1080 c call hrget(nlamx,'hbooks/ngcx.hbook','') c thx0 = HRNDM1(nlamx) c call hdelet(nlamx) c call hrget(nlamy,'hbooks/ngcy.hbook','') c thy0 = HRNDM1(nlamy) c call hdelet(nlamy) C************************************************************** CC**************************************************************** c thxinput = thx0 thyinput = thy0 call hf1(31,real(thxinput),1.) call hf1(32,real(thyinput),1.) NPSMout = NPSMout + 1 ! bec count exiting PSM PSMAVGflux = PSMAVGflux + 1/xin/yin/4.0d0/nevts if (x0**2+y0**2.le.0.625**2) & PSMPKflux = PSMPKflux + 1/pi/0.625**2/nevts capflpsm = capflpsm + sqrt(25.0)/rsqrt(en0) call hf1(1,real(x0),1.) call hf1(10,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) ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp 2341 format(' EastDn ',1x,f9.2,1x,f9.2,1x,f9.2) 2342 format(' WestUp ',1x,f9.2,1x,f9.2,1x,f9.2) ******************************************************************** * Neutrons passing through the gap (3.0 cm) after the PSM before input guides * this is being added in to old code, so x0, y0, z0, etc just redefined as * result of this 3 cm transport * if((x0.ge.tairx1.and.x0.le.tairx2).and. & (y0.ge.tairy1.and.y0.le.tairy2)) Nairin1 = Nairin1 + 1 c call airgap(en0, x0, y0, z0, thx0, thy0, & tairx1,tairx2,tairy1,tairy2,tairz1, en1, x1, y1, z1, & thx1,thy1,airflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (airflg1) Nairscat1 = Nairscat1+1 ! count scatters with energy change bec if (airflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(airflg1)).or. & ((en0.ne.en1).and.(.not.(airflg1))) ) & write(6,*) ' Edif/airflg AIR 1',airflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif AIR 1 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg c write(6,*) x1,airflg1,beamscflg if (en1.lt.0.d0) then write(6,*) 'ag1 en<1, en=1e-6' en1 = 0.0000001 end if pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot ! mrad if (pltot.gt.1300) write(6,*) '0th air gap',pltot Nairout1 = Nairout1 + 1 if (dabs(thtot).gt.1.0e9) then write(6,*) 'Air thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** * Neutrons passing through FIRST Al Window * 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 if((x0.ge.talx1.and.x0.le.talx2).and. & (y0.ge.taly1.and.y0.le.taly2)) Nalin1 = Nalin1 + 1 c call alwindow(en0, x0, y0, z0, thx0, thy0, & talx1,talx2,taly1,taly2,talz1, en1, x1, y1, z1, & thx1,thy1,alflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (alflg1) Nalscat1 = Nalscat1+1 ! count scatters with energy change bec if (alflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(alflg1)).or. & ((en0.ne.en1).and.(.not.(alflg1))) ) & write(6,*) ' Edif/airflg AL 1',alflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif AIR 1 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg if (en1.lt.0.d0) then write(6,*) 'Al1 en<1, en=1e-6' en1 = 0.0000001 end if pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.1300) write(6,*) 'Al 1 ',pltot Nalout1 = Nalout1 + 1 if (dabs(thtot).gt.1.0e9) then write(6,*) ' Al thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** c Neutrons passing through the guide * 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 beamflg = 1 if (sptmflg) then if ((x0.ge.gx1.and.x0.le.gx2) & .and.(y0.ge.gy1.and.y0.lt.gy2)) then call beampipe(beamflg,en0,x0,y0,z0, & thx0,thy0,gx1,gx2,gy1,gy2,zlm,en1,x1,y1,z1, & thx1,thy1,scflag,nbx,nby,pl_del,rotang) if((x0*x1.lt.0d0).and.x1.gt.-1d5) & write(6,*) ' Guide Crossover ',x0,x1 !check for crossovers else if ((x0.ge.gx3.and.x0.le.gx4).and. & (y0.ge.gy1.and.y0.lt.gy2)) then call beampipe(beamflg,en0,x0,y0,z0, & thx0,thy0,gx3,gx4,gy1,gy2,zlm,en1,x1,y1,z1, & thx1,thy1,scflag,nbx,nby,pl_del,rotang) if((x0*x1.lt.0d0).and.x1.gt.-1d5) & write(6,*) ' Guide Crossover ',x0,x1 !check for crossovers else goto 100 end if else if((x0.ge.gx1.and.x0.le.gx4).and.(y0.ge.gy1.and.y0.le.gy2)) then call beampipe(beamflg,en0,x0,y0,z0, & thx0,thy0,gx1,gx4,gy1,gy2,zlm,en1,x1,y1,z1, & thx1,thy1,scflag,nbx,nby,pl_del,rotang) else goto 100 end if end if if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 call hf1(2,real(x1),1.) if (nbx.eq.0) ngbx0=ngbx0+1 if (nbx.eq.1) ngbx1=ngbx1+1 if (nbx.eq.2) ngbx2=ngbx2+1 if (nbx.eq.3) ngbx3=ngbx3+1 if (nby.eq.0) ngby0=ngby0+1 if (nby.eq.1) ngby1=ngby1+1 if (nby.eq.2) ngby2=ngby2+1 if (nby.eq.3) ngby3=ngby3+1 if (nbx.gt.3) write(6,*) ' 4+ Guide xBounces ',nbx if (nby.gt.3) write(6,*) ' 4+ Guide yBounces ',nby gsc = gsc + scflag if(scflag.gt.5) write(6,*)' 5+ Guide Scatters ',scflag if (scflag.gt.0) beamscflg = .true. if ( ((en0.eq.en1).and.(scflag.gt.0)).or. & ((en0.ne.en1).and.(scflag.eq.0)) ) & write(6,*) ' Edif/scatflg GUIDE',scflag,en0,en1 c if (en0.ne.en1) write(6,*) ' edif AIR 1 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg c This path length is now OK for reflected neutrons -KG c pl_del = gl/dcos(ratan(rsqrt((dtan(thx0))**2+(dtan(thy0))**2))) c pltot = pl_del + pltot c thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot pltot = pl_del + pltot thtot = thtot + rotang if (pltot.gt.1300) write(6,*) 'guide ',pltot Nguideout = Nguideout+1 if(x1.gt.gx1.and.x1.lt.gx2) NguideoutE = NguideoutE+wgt if(x1.gt.gx3.and.x1.lt.gx4) NguideoutW = NguideoutW+wgt if (dabs(thtot).gt.1.0e9) then write(6,*) ' guide thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if if ((thx1.eq.-1.0*thx0).or.(thy1.eq.-1.0*thy0)) bounceflg = .true. ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** ******************************************************************** * Neutrons passing through SECOND Al Window * 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 if((x0.ge.talx1.and.x0.le.talx2).and. & (y0.ge.taly1.and.y0.le.taly2)) Nalin2 = Nalin2 + 1 ! rcm c call alwindow(en0, x0, y0, z0, thx0, thy0, & talx1,talx2,taly1,taly2,talz2, en1, x1, y1, z1, & thx1,thy1,alflg2) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (alflg2) Nalscat2 = Nalscat2+1 !rcm if (alflg2) beamscflg = .true. if ( ((en0.eq.en1).and.(alflg1)).or. & ((en0.ne.en1).and.(.not.(alflg1))) ) & write(6,*) ' Edif/alflg Al 2',alflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif AIR 1 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg if(en1.lt.0.d0) then write(6,*) 'All en<1, en=1e-6' en1 = 0.0000001 end if pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.1300) write(6,*) 'Al 2 ',pltot Nalout2 = Nalout2 + 1 if (dabs(thtot).gt.1.0e9) then write(6,*) ' Al thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ************************************************************************* * Neutrons passing through Al wire input coil entrance * 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 alwindow(en0, x0, y0, z0, thx0, thy0, & talx1,talx2,taly1,taly2,talz6, en1, x1, y1, z1, & thx1,thy1,alflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (alflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(alflg1)).or. & ((en0.ne.en1).and.(.not.(alflg1))) ) & write(6,*) ' Edif/alflg Al 3',alflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif AIR 1 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.1300) write(6,*) 'Al input ',pltot if (dabs(thtot).gt.1.0e9) then write(6,*) ' Al out thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ************************************************************************ * Neutrons passing through the input coil * 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 beamflg = 2 if (sptmflg) then if ((x0.ge.ix1.and.x0.le.ix2).and. & (y0.ge.iy1.and.y0.le.iy2)) then call beampipe(beamflg,en0,x0,y0,z0, & thx0,thy0,ix1,ix2,iy1,iy2,zlm,en1,x1,y1,z1, & thx1,thy1,scflag,nbx,nby,pl_del,rotang) if((x0*x1.lt.0d0).and.x1.gt.-1d5) & write(6,*) ' Incoil Crossover ',x0,x1 !check for crossovers else if ((x0.ge.ix3.and.x0.le.ix4) & .and.(y0.ge.iy1.and.y0.le.iy2)) then call beampipe(beamflg,en0,x0,y0,z0, & thx0,thy0,ix3,ix4,iy1,iy2,zlm,en1,x1,y1,z1, & thx1,thy1,scflag,nbx,nby,pl_del,rotang) if((x0*x1.lt.0d0).and.x1.gt.-1d5) & write(6,*) ' Incoil Crossover ',x0,x1 !check for crossovers else goto 100 end if else if((x0.ge.ix1.and.x0.le.ix4).and.(y0.ge.iy1.and.y0.le.iy2)) then call beampipe(beamflg,en0,x0,y0,z0, & thx0,thy0,ix1,ix4,iy1,iy2,zlm,en1,x1,y1,z1, & thx1,thy1,scflag,nbx,nby,pl_del,rotang) else goto 100 end if end if if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 Nincoilout = Nincoilout + 1 ! bec counttransmission through guide and Input coil if(x1.gt.ix1.and.x1.lt.ix2) NincoiloutE = NincoiloutE+wgt if(x1.gt.ix3.and.x1.lt.ix4) NincoiloutW = NincoiloutW+wgt call hf1(3,real(x1),1.) call hf2(1111,real(x1),real(y1),1.) if (nbx.eq.0) nibx0=nibx0+1 if (nbx.eq.1) nibx1=nibx1+1 if (nbx.eq.2) nibx2=nibx2+1 if (nbx.eq.3) nibx3=nibx3+1 if (nby.eq.0) niby0=niby0+1 if (nby.eq.1) niby1=niby1+1 if (nby.eq.2) niby2=niby2+1 if (nby.eq.3) niby3=niby3+1 if (nbx.gt.3) write(6,*) ' 4+ Incoil xBounces ',nbx if (nby.gt.3) write(6,*) ' 4+ Incoil yBounces ',nby isc = isc + scflag if(scflag.gt.5) write(6,*)' 5+ Incoil Scatters ',scflag if (scflag.gt.0) beamscflg = .true. if ( ((en0.eq.en1).and.(scflag.gt.0)).or. & ((en0.ne.en1).and.(scflag.eq.0)) ) & write(6,*) ' Edif/scatflg INPUT',scflag,en0,en1 c if (en0.ne.en1) write(6,*) ' edif INCOIL ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg c This path length is now OK for reflected neutrons -KG cc pl_del = iz0/dcos(ratan(rsqrt((dtan(thx0))**2+(dtan(thy0))**2))) bec2011 no longer needed with new routines! pltot = pl_del + pltot thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.1300) write(6,*) 'incoil ',pltot if (dabs(thtot).gt.1.0e9) then write(6,*) ' incoil thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if * if ((bounceflg).and.((thx1.eq.-1.0*thx0).or. & (thy1.eq.-1.0*thy0))) bounceflg2 = .true. if ((thx1.eq.-1.0*thx0).or.(thy1.eq.-1.0*thy0)) bounceflg = .true. capflinc = capflinc + sqrt(25.0)/rsqrt(en1) ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** * Neutrons passing through Al wire input coil exit * 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 alwindow(en0, x0, y0, z0, thx0, thy0, & talx1,talx2,taly1,taly2,talz6, en1, x1, y1, z1, & thx1,thy1,alflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (alflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(alflg1)).or. & ((en0.ne.en1).and.(.not.(alflg1))) ) & write(6,*) ' Edif/alflg Al 3',alflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif Al inexit ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.1300) write(6,*) 'Al input ',pltot if (dabs(thtot).gt.1.0e9) then write(6,*) ' Al out thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** ******************************************************************** * 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 (airflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(airflg1)).or. & ((en0.ne.en1).and.(.not.(airflg1))) ) & write(6,*) ' Edif/airflg Al 3',airflg1,en0,en1 c call collimator(x0,x1,y0,y1) !Boroflex septum to prevent crossovers c if (en0.ne.en1) write(6,*) ' edif in/tar ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg 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_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.1300) 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 * ******************************************************************** ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************* INput Coil Exit Flux Measurement*************** if (sptmflg) then INCOAVGflux = INCOAVGflux + 1/((ix2-ix1)+(ix4-ix3))/(iy2-iy1)/nevts if ((x1+(ix2+ix1)/2.)**2+y1**2.le.0.625**2) & INCOPKflux = INCOPKflux + 1/pi/0.625**2/nevts else INCOAVGflux = INCOAVGflux + 1/(ix4-ix1)/(iy2-iy1)/nevts if (x1**2+y1**2.le.0.625**2) & INCOPKflux = INCOPKflux + 1/pi/0.625**2/nevts end if 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) ******************************************************************** ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** * 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 if (airflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(airflg1)).or. & ((en0.ne.en1).and.(.not.(airflg1))) ) & write(6,*) ' Edif/airflg AIR in/tar',airflg1,en0,en1 c call collimator(x0,x1,y0,y1) !Boroflex septum to prevent crossovers c if (en0.ne.en1) write(6,*) ' edif in/tar2 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg 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_beam*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.1300) 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 * ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** * Neutrons passing through THIRD Al window * 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 alwindow(en0, x0, y0, z0, thx0, thy0, & talx1,talx2,taly1,taly2,talz3, en1, x1, y1, z1, & thx1,thy1,alflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (alflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(alflg1)).or. & ((en0.ne.en1).and.(.not.(alflg1))) ) & write(6,*) ' Edif/alflg Al 3',alflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif Al 3 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg if (en1.lt.0.d0) then write(6,*) 'Al2 en<1, en=1e-6' en1 = 0.0000001 end if if (alflg1) beamscflg = .true. pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.1300) write(6,*) 'Al 3 ',pltot if (dabs(thtot).gt.1.0e9) then write(6,*) ' Al thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if ******************************************************************** NtarWEin = NtarWEin + 1 xtar1 = x1 ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** c print *, " first half target" * Neutrons passing through the 1st half of the target (41.6 cm). C Note: below checks only where the neutron is in x direction. C Once it enters target target.f checks the y-direction too, so might C have neutrons that get into target.f that don't run target.f really, C since the y-values are out of range, hence sig or thscat values get C nothing put in them * * Defining x<0 to imply West which in main_WUp.exe has upstream target C c Update x,y,z,thx,thy,en values to new values x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1 ztar1=z0 c Given current energy determine q from S(q) and delE from omega(q) k0 = rsqrt(2*m_n*en0)/h_par c c enflg1 = .false. thflg1 = .false. empflg1 = 0 if( (x0.ge.tarx1.and.x0.le.tarx2) ) then ! East=neg. x, empty UPstream if (y0.ge.tary1.and.y0.le.tary2) Nemp1in=Nemp1in+wgt call empty(en0, x0, y0, z0, thx0, thy0, tarx1, tarx2, & en1, x1, y1, z1, thx1, thy1,empflg1, ztar1,pl_del,rotang) if (en1.lt.0.d0) goto 100 c write(6,*) 'tar1 en<1, en=1e-6' c en1 = 0.0000001 c end if c write(6,*) ' 1st emp E ',z1,pl_del if (empflg1.gt.0) Nempty1 = Nempty1+1.0 ! not trajectory dependent -- no wgt if (empflg1.gt.0) beamscflg = .true. if ((x1.ge.tarx1.and.x1.le.tarx2).and. & (y1.ge.tary1.and.y1.le.tary2)) & Nemp1out = Nemp1out+wgt ! out of empty west side only elseif( (x0.ge.tarx3.and.x0.le.tarx4) ) then ! West=pos. x, Target UP sig = 0.0 ; thscat = 0.0 thp = 0.0; php = 0.0 loop = 0 qmin = 0 qmax = int(k0*1e3) if(qmax.lt.100) then hmax = qmax else hmax=100 endif call hcopyr(699,302,"",qmin,hmax,0.,0.,"") cbec10 call hcopyr(700,302,"",qmin,hmax,0.,0.,"") 501 q=HRNDM1(302) if(q.eq.0.) goto 501 hws1 = q*1e3+700 call hrget(hws1,'hbooks/sqw.hbook','') wp1=HRNDM1(hws1) call hdelet(hws1) ep1 = en0 - (h_par*wp1) if (ep1.lt.0.d0) then write(6,*) ' 1st tar, ep1<0 ',ep1 goto 501 endif kp1 = rsqrt(2*m_n*ep1)/h_par th_test1 = (k0**2 + kp1**2 - q**2)/(2.*k0*kp1) loop = loop + 1 if(ep1.lt.0.or.abs(th_test1).gt.1.0) goto 501 call hdelet(302) call hf1(305,real(q),1.) call hf1(306,real(wp1),1.) call hf1(307,real(en0),1.) call hf1(308,real(ep1),1.) call hf1(309,real(k0),1.) call hf1(310,real(kp1),1.) call target(q, en0, x0, y0, z0, thx0, thy0, tarx3, tarx4, & en1, x1, y1, z1, thx1, thy1, sig, thscat, & scatflg1,nrfl1,ztar1,pl_del,rotang,scwgt) ! bec if (en1.lt.-1000.d0) goto 100 c make sure really went in target -- we know x-position is ok, check y-position if (y0.ge.tary1.and.y0.le.tary2) then c write(6,*) ' 1st tar W ',z1,pl_del if(en1.ne.en0) enflg1=.true. If ((.not.hegastar).and.(.not.airtar)) wgt2=wgtA/sig*scwgt c If ((.not.hegastar).and.(.not.airtar)) wgt2=1.0d0 c if ((en1.gt.-1d5).and.(abs(wgt2).gt.1)) c & write(*,*) ' wgt2 ',wgt2, en0, en1, sig ! bec 2010 Ntar1=Ntar1+wgt ! entering target, before wgt change cbec fall2011 if((en0.gt.-100).and.(en1.gt.-100).and.(en0.ne.en1)) then ! scattered if (scatflg1.gt.0) then ! scattered in target.f/scatt.f tar1flg= .true. c write(81,*) en0,en1,en0-en1 c determine the scattering weighting factor ratio of integral of S(q,w) c over limited q range (q<0.1 1/Ang) to full c integral (= energy dependent cross section from Sommers) If ((.not.hegastar).and.(.not.airtar)) wgt=wgtA/sig*scwgt c If ((.not.hegastar).and.(.not.airtar)) wgt=1.0d0 call hf1(7001,real(en0-en1),wgt) endif if ( (enflg1.and.(scatflg1.eq.0)).or. & (.not.(enflg1).and.(scatflg1.gt.0)) ) & write(6,*) ' enflg/scatflg ',scatflg1,en0,en1 if (wgt.lt.0) write(6,*) ' 1st tar wgt<0 ', wgt, sig, en1 if (enflg1) Nenscat1 = Nenscat1+1.0 ! not trajectory dependent -- no wgt c if (thflg1) Nthscat1 = Nthscat1+1.0 ! not trajectory dependent -- no wgt c capture diagnostics for target interaction stuff bec call hf1(4010,real(sig),1.0) !not trajectory dependent -- no wgt call hf2(8010,real(q),real(sig),1.0) !temporary!not trajectory dependent -- no wgt call hf1(4012,real(180./3.14159*thscat),wgt) call hf1(4015,real(180./3.14159*thp),wgt) c call hf1(4016,real(180./3.14159*php),wgt) end if if ((x1.ge.tarx3.and.x1.le.tarx4).and. & (y1.ge.tary1.and.y1.le.tary2)) & Ntar1out = Ntar1out+wgt ! out of target west side only else x1 = -1.d6 y1 = -1.d6 z1 = -1.d6 thx1 = -1.d6 thy1 = -1.d6 en1 = -1.d6 endif * if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 * entar1 = en1 call hf1(1001,real(x1),wgt) call hf1(1002,real(y1),wgt) c call hf1(1003,real(thx1),wgt) c call hf1(1004,real(thy1),wgt) cccccccccccccc target rotation c write(6,*) ' pl_del ', pl_del, ' rotang ',rotang pltot = pl_del + pltot thtot = thtot + rotang cbec fall2011 if(en0.eq.en1) then if (scatflg1.eq.0) then call hf1(1005,real(en1),wgt) else call hf1(1012,real(en1),wgt) endif cccccccccc end target rotation if (dabs(thtot).gt.1.0e9) then write(6,*) ' tar1 thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if c write(*,*) 'tar 1 ends at z = ',z1 !-KG * ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** * Neutrons passing through FOURTH Al window * 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 alwindow(en0, x0, y0, z0, thx0, thy0, & talx1,talx2,taly1,taly2,talz4, en1, x1, y1, z1, & thx1,thy1,alflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (alflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(alflg1)).or. & ((en0.ne.en1).and.(.not.(alflg1))) ) & write(6,*) ' Edif/alflg Al 4',alflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif Al 4 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg call rotation(en0,x0,y0,z0,x1,y1,z1,ztar1,rotang) ! ztar1 is tarting point of gradient field if being used pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = thtot + rotang if (pltot.gt.1300) write(6,*) 'Al 4 ',pltot if (dabs(thtot).gt.1.0e9) then write(6,*) ' Al tar1 thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if ******************************************************************** ctemp if (wgt.ne.1.0) write(17,*) ' tar1 out wgt ',wgt,tar1flg,tar2flg ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** * Neutrons passing through Al pi-coil wire support * 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 alwindow(en0, x0, y0, z0, thx0, thy0, & talx1,talx2,taly1,taly2,talz5, en1, x1, y1, z1, & thx1,thy1,alflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (alflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(alflg1)).or. & ((en0.ne.en1).and.(.not.(alflg1))) ) & write(6,*) ' Edif/alflg Al pi1',alflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif Al pi1 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg call rotation(en0,x0,y0,z0,x1,y1,z1,ztar1,rotang) ! ztar1 is tarting point of gradient field if being used pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = thtot + rotang if (pltot.gt.1300) write(6,*) 'Al Pi ',pltot,x0,y0,z0,x1,y1,z1 if (dabs(thtot).gt.1.0e9) then write(6,*) ' Al out thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp *********************************************************************** c print *," Gap btw targets" * Neutrons passing thruogh the gap (pi-coil = 7.3 cm) between the * two tagets = 8.0cm. * 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 cccc 1st half of air gap call airgap(en0, x0, y0, z0, thx0, thy0, & tairx1,tairx2,tairy1,tairy2,0.5*tairz3, en1, x1, y1, z1, & thx1,thy1,airflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (airflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(airflg1)).or. & ((en0.ne.en1).and.(.not.(airflg1))) ) & write(6,*) ' Edif/airflg AIR tar1',airflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif AIR tar1 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg call rotation(en0,x0,y0,z0,x1,y1,z1,ztar1,rotang) ! ztar1 is starting point of gradient field if being used pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = thtot + rotang call hf1(5019,real(180./3.14159*thx1),wgt) call hf1(5020,real(180./3.14159*thy1),wgt) 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(5022,real(180./3.14159*phi),wgt) ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ********************************************************************* * Neutrons passing through Al Wire Pi Coil support * 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 alwindow(en0, x0, y0, z0, thx0, thy0, & talx1,talx2,taly1,taly2,talz5, en1, x1, y1, z1, & thx1,thy1,alflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (alflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(alflg1)).or. & ((en0.ne.en1).and.(.not.(alflg1))) ) & write(6,*) ' Edif/alflg Al pi2',alflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif Al pi2 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg call rotation(en0,x0,y0,z0,x1,y1,z1,ztar1,rotang) ! ztar1 is tarting point of gradient field if being used pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = thtot + rotang if (pltot.gt.1300) write(6,*) 'Al Pi ',pltot,x0,y0,z0,x1,y1,z1 if (dabs(thtot).gt.1.0e9) then write(6,*) ' Al out thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp *********************************************************************** c do NOT allow neutrons from one target cell into opposite side cell if( ((x0.ge.tarx1).and.(x0.le.tarx2)).and.((x1.lt.tarx1).or. & (x1.gt.tarx2)) ) & goto 100 if( ((x0.ge.tarx3).and.(x0.le.tarx4)).and.((x1.lt.tarx3).or. & (x1.gt.tarx4)) ) & goto 100 c*************PI Coil reversal of rotation angle c using 1/2 gap length -- rotate to halfway point, reverse, rotate the other half c reverse thtot in center of gap and rotate the other half of gap length if (picoilflg) then ! pi-coil ON thtot = - thtot end if c Update x,y,z,thx,thy,en values to new values x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1 call airgap(en0, x0, y0, z0, thx0, thy0, & tairx1,tairx2,tairy1,tairy2,0.5*tairz3, en1, x1, y1, z1, & thx1,thy1,airflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (airflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(airflg1)).or. & ((en0.ne.en1).and.(.not.(airflg1))) ) & write(6,*) ' Edif/airflg AIR tar2',airflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif AIR tar2 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg call rotation(en0,x0,y0,z0,x1,y1,z1,ztar1,rotang) ! ztar1 is tarting point of gradient field if being used pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = thtot + rotang if (dabs(thtot).gt.1.0e9) then write(6,*) ' thtot ', thtot, ', pltot ',pltot write(6,*) ' en0 ',en0,', en ',en1 write(6,*) ' thxy01',thx0,thy0, thx1, thy1 end if ctemp cc if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** * Neutrons passing through FIFTH Al window * 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 alwindow(en0, x0, y0, z0, thx0, thy0, & talx1,talx2,taly1,taly2,talz4, en1, x1, y1, z1, & thx1,thy1,alflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (alflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(alflg1)).or. & ((en0.ne.en1).and.(.not.(alflg1))) ) & write(6,*) ' Edif/alflg Al 5',alflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif Al 5 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg call rotation(en0,x0,y0,z0,x1,y1,z1,ztar1,rotang) ! ztar1 is tarting point of gradient field if being used pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = thtot + rotang if (pltot.gt.1300) write(6,*) 'Al 5 ',pltot if (dabs(thtot).gt.1.0e9) then write(6,*) ' Al thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if ******************************************************************** xtar2 = x1 cctemp if (wgt.ne.1.0) write(17,*) ' tar2 in wgt ',wgt,tar1flg,tar2flg ctemp cc if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** c print *, " second half target" * Neutrons passing thruogh the 2nd part of the target (41.6 cm). * * 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 enflg2 = .false. thflg2 = .false. empflg2 = 0 if( (x0.gt.tarx1.and.x0.lt.tarx2) ) then ! East = neg. x, Target down if (tar1flg) then ! if scat in t1, now entered both targets Ntar12=Ntar12+wgt endif sig = 0.0 ; thscat = 0.0 thp = 0.0; php = 0.0 c c Given current energy determine q from S(q) and delE from omega(q) c Doing this only if in east side where target is located c k0 = rsqrt(2*m_n*en0)/h_par loop = 0 qmin = 0 qmax = int(k0*1e3) if(qmax.lt.100) then hmax = qmax else hmax=100 endif call hcopyr(699,302,"",qmin,hmax,0.,0.,"") cbec10 call hcopyr(700,302,"",qmin,hmax,0.,0.,"") 502 q=HRNDM1(302) if(q.eq.0.) goto 502 hws2 = q*1e3+700 call hrget(hws2,'hbooks/sqw.hbook','') wp2=HRNDM1(hws2) call hdelet(hws2) ep2 = en0 - (h_par*wp2) if (ep2.lt.0.d0) then write(6,*) ' 2nd tar, ep2<0 ',ep2 goto 502 endif kp2 = rsqrt(2*m_n*ep2)/h_par th_test2 = (k0**2 + kp2**2 - q**2)/(2.*k0*kp2) loop = loop + 1 if(ep2.lt.0.or.abs(th_test2).gt.1.0) goto 502 call hdelet(302) call hf1(205,real(q),wgt) call hf1(206,real(wp2),wgt) call hf1(207,real(en0),wgt) call hf1(208,real(ep2),wgt) call hf1(209,real(k0),wgt) call hf1(210,real(kp2),wgt) call target(q, en0, x0, y0, z0, thx0, thy0, tarx1, tarx2, & en1, x1, y1, z1, thx1, thy1, sig, thscat, & scatflg2, nrfl2, ztar1, pl_del, rotang,scwgt) ! bec if (en1.lt.-1000.d0) goto 100 c make sure really went in target -- we know x-position is ok, check y-position if (y0.ge.tary1.and.y0.le.tary2) then if (en1.ne.en0) enflg2 = .true. If ((.not.hegastar).and.(.not.airtar)) wgt2=wgtA/sig*scwgt c If ((.not.hegastar).and.(.not.airtar)) wgt2=1.0d0 Ntar2=Ntar2+wgt ! entering target, before wgt change c cbec fall2011 if(en0.gt.-100.and.en1.gt.-100.and.en0.ne.en1) then ! scattered if (scatflg2.gt.0) then tar2flg = .true. c write(81,*) en0,en1,en0-en1 c determine the scattering weighting factor ratio of integral of S(q,w) c over limited q range (q<0.1 1/Ang) to full c integral (= energy dependent cross section from Sommers) If ((.not.hegastar).and.(.not.airtar)) wgt=wgtA/sig*scwgt c If ((.not.hegastar).and.(.not.airtar)) wgt=1.0d0 call hf1(7002,real(en0-en1),wgt) endif if ( (enflg2.and.(scatflg2.eq.0)).or. & (.not.(enflg2).and.(scatflg2.gt.0)) ) & write(6,*) ' enflg/scatflg ',scatflg2,en0,en1 if (wgt.lt.0) write(6,*) ' 2nd tar wgt<0 ',wgt, sig, en1 c if (enflg2) Nenscat2 = Nenscat2+1.0 ! not trajectory dependent -- no wgt c if (thflg2) Nthscat2 = Nthscat2+1.0 ! not trajectory dependent -- no wgt c capture diagnostics for target interaction stuff bec call hf1(4020,real(sig),1.0) c call hf1(4021,real(intlen),wgt) call hf1(5012,real(180./3.14159*thscat),wgt) call hf1(5015,real(180./3.14159*thp),wgt) c call hf1(5016,real(180./3.14159*php),wgt) call hf2(9998,real(q),real(wp1),wgt) ! Dispersion relation -KG end if if ((x1.ge.tarx1.and.x1.le.tarx2).and. & (y1.ge.tary1.and.y1.le.tary2)) Ntar2out = Ntar2out+wgt ! out of target west side only elseif( (x0.gt.tarx3.and.x0.lt.tarx4) ) then ! West = pos. x, empty DOWN if (y0.ge.tary1.and.y0.le.tary2) Nemp2in=Nemp2in+wgt call empty(en1, x0, y0, z0, thx0, thy0, tarx3, tarx4, & en1, x1, y1, z1, thx1, thy1, empflg2,ztar1,pl_del,rotang) if (en1.lt.0.d0) goto 100 if (empflg2.gt.0) Nempty2 = Nempty2+1.0 ! not trajectory dependent -- no wgt if (empflg2.gt.0) beamscflg = .true. if ((x1.ge.tarx3.and.x1.le.tarx4).and. & (y1.ge.tary1.and.y1.le.tary2)) & Nemp2out = Nemp2out+wgt ! out of empty east side only else x1 = -1.d6 y1 = -1.d6 z1 = -1.d6 thx1 = -1.d6 thy1 = -1.d6 en1 = -1.d6 endif * if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 * entar2 = en1 call hf1(2001,real(x1),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) c call hf1(2003,real(thx1),wgt) c call hf1(2004,real(thy1),wgt) cccccccc target 2 rotation pltot = pl_del + pltot thtot = thtot + rotang cbec fall2011 if(en1.eq.en0) then if (scatflg2.eq.0) then call hf1(2005,real(en1),wgt) else call hf1(2012,real(en1),wgt) endif cccccccccccc end target 2 rotaton if (pltot.gt.1300) write(6,*) 'target 2 ',pltot if (dabs(thtot).gt.1.0e9) then write(6,*) ' tar2 thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if c write(*,*) 'tar 2 ends at z = ', z1 * ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot cc if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** ctemp if (wgt.ne.1.0) write(17,*) ' tar2 out wgt ',wgt,tar1flg,tar2flg NTargetout = NTargetout + wgt TARAVGflux = TARAVGflux + & wgt/((tarx2-tarx1)+(tarx4-tarx3))/(tary2-tary1)/nevts if ((x1-0.5*(tarx2-tarx1))**2+y1**2.le.0.625**2) & TARPKflux = TARPKflux + wgt/pi/0.625**2/nevts ******************************************************************** ******************************************************************** * Neutrons passing through SIXTH Al window * 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 alwindow(en0, x0, y0, z0, thx0, thy0, & talx1,talx2,taly1,taly2,talz3, en1, x1, y1, z1, & thx1,thy1,alflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (alflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(alflg1)).or. & ((en0.ne.en1).and.(.not.(alflg1))) ) & write(6,*) ' Edif/alflg Al 6',alflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif Al 6 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg 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.1300) write(6,*) 'Al 6 ',pltot if (dabs(thtot).gt.1.0e9) then write(6,*) ' Al tar2 thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if ctemp cc if (x1.lt.0) write(6,2341) z1,pltot,thtot cc if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp *********************************************************************** * Neutrons passing thruogh the gap (13.0 cm) between the taget and * the output coil. * 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,tairz4, 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 if (airflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(airflg1)).or. & ((en0.ne.en1).and.(.not.(airflg1))) ) & write(6,*) ' Edif/airflg AIR tar/out',airflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif AIR tar/out ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg 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(5023,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_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.1300) write(6,*) 'gap 3 ',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 ******************************************************************** ctemp cc if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** * Neutrons passing through Al Wire output coil entrance * 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 alwindow(en0, x0, y0, z0, thx0, thy0, & talx1,talx2,taly1,taly2,talz6, en1, x1, y1, z1, & thx1,thy1,alflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (alflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(alflg1)).or. & ((en0.ne.en1).and.(.not.(alflg1))) ) & write(6,*) ' Edif/alflg Al outcoil in',alflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif Al out1 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg 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.1300) write(6,*) 'Al output ',pltot if (dabs(thtot).gt.1.0e9) then write(6,*) ' Al out thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** * Neutrons passing through the output coil. * c Update x,y,z,thx,thy,en values to new values x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1 cccccccccccccc If you want to put outcoil guide at an angle ccccc c if (i.eq.1) write(6,*) ' !!! outcoil ang offset, thy->2mrad ' c thy0 = thy0+0.002 cccccccccccccccccccccc beamflg = 3 if((x0.ge.ox1.and.x0.le.ox4).and.(y0.ge.oy1.and.y0.le.oy2)) then if (x0.ge.ox1.and.x0.le.ox2) then call beampipe(beamflg,en0,x0,y0,z0, & thx0,thy0,ox1,ox2,oy1,oy2,zlm,en1,x1,y1,z1, & thx1,thy1,scflag,nbx,nby,pl_del,rotang) if((x0*x1.lt.0d0).and.x1.gt.-1d5) & write(6,*) ' Outcoil Crossover ',x0,x1 !check for crossovers else if (x0.ge.ox3.and.x0.le.ox4) then call beampipe(beamflg,en0,x0,y0,z0, & thx0,thy0,ox3,ox4,oy1,oy2,zlm,en1,x1,y1,z1, & thx1,thy1,scflag,nbx,nby,pl_del,rotang) if((x0*x1.lt.0d0).and.(x1.gt.-1d5)) & write(6,*) 'Outcoil Crossover ',x0,x1 !check for crossovers else goto 100 end if else goto 100 end if 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 (nbx.eq.0) nobx0=nobx0+1 if (nbx.eq.1) nobx1=nobx1+1 if (nbx.eq.2) nobx2=nobx2+1 if (nbx.eq.3) nobx3=nobx3+1 if (nby.eq.0) noby0=noby0+1 if (nby.eq.1) noby1=noby1+1 if (nby.eq.2) noby2=noby2+1 if (nby.eq.3) noby3=noby3+1 if (nbx.gt.3) write(6,*) ' 4+ Outcoil xBounces ',nbx if (nby.gt.3) write(6,*) ' 4+ Outcoil yBounces ',nby osc = osc + scflag if(scflag.gt.5) write(6,*)' 5+ Outcoil Scatters ',scflag if (scflag.gt.0) beamscflg = .true. if ( ((en0.eq.en1).and.(scflag.gt.0)).or. & ((en0.ne.en1).and.(scflag.eq.0)) ) & write(6,*) ' Edif/scatflg INPUT',scflag,en0,en1 c if (en0.ne.en1) write(6,*) ' edif out ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg c This path length is not OK for reflected neutrons -KG bec2011 no longer needed with new routines c pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) c pl_del = oz/dcos(ratan(rsqrt((dtan(thx0))**2+(dtan(thy0))**2))) pltot = pl_del + pltot thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.1300) write(6,*) 'out coil ',pltot if (dabs(thtot).gt.1.0e9) then write(6,*) ' thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot end if c if ((bounceflg).and.(bounceflg2).and.((thx1.eq.-1.0*thx0).or. & (thy1.eq.-1.0*thy0))) bounceflg3 = .true. if ((bounceflg).and.((thx1.eq.-1.0*thx0).or. & (thy1.eq.-1.0*thy0))) bounceflg2 = .true. if ((thx1.eq.-1.0*thx0).or.(thy1.eq.-1.0*thy0)) bounceflg = .true. capflout = capflout + sqrt(25.0)/rsqrt(en1) * ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot cc if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** * Neutrons passing through Al Wire output coil exit * 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 alwindow(en0, x0, y0, z0, thx0, thy0, & talx1,talx2,taly1,taly2,talz6, en1, x1, y1, z1, & thx1,thy1,alflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (alflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(alflg1)).or. & ((en0.ne.en1).and.(.not.(alflg1))) ) & write(6,*) ' Edif/alflg Al outcoil',alflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif Al out2 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.1300) write(6,*) 'Al output ',pltot if (dabs(thtot).gt.1.0e9) then write(6,*) ' Al out thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if ************************************************************************ ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ************************************************************************ * Neutrons passing through the gap (10.0 cm) between the output coil * and the analyzing super mirror (ASM). * 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,tairz5, 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 if (airflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(airflg1)).or. & ((en0.ne.en1).and.(.not.(airflg1))) ) & write(6,*) ' Edif/airflg Al 3',airflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif AIR out/asm ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg NOutcoilout = NOutcoilout + wgt if(x1.gt.ox1.and.x1.lt.ox2) NoutcoiloutE = NoutcoiloutE+wgt if(x1.gt.ox3.and.x1.lt.ox4) NoutcoiloutW = NoutcoiloutW+wgt OUTCOAVGflux = OUTCOAVGflux + & wgt/((ox2-ox1)+(ox4-ox3))/(oy2-oy1)/nevts if (x1**2+y1**2.le.0.625**2) & OUTCOPKflux = OUTCOPKflux + wgt/pi/0.625**2/nevts pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.1300) write(6,*) 'gap 5 ',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 *************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) ******************************************************************** ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** * Neutrons passing through the super mirror." c print *," Super Mirror" * 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 sm(en0,x0,y0,z0,thx0,thy0,x1,y1,z1,thx1,thy1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 cot=cot+wgt call hf1(6,real(x1),1.) pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.1300) write(6,*) 'ASM ',pltot if (dabs(thtot).gt.1.0e9) then write(6,*) ' ASM thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if * if ((bounceflg).and.(bounceflg2).and.((thx1.eq.-1.0*thx0).or. & (thy1.eq.-1.0*thy0))) bounceflg3 = .true. if ((bounceflg).and.((thx1.eq.-1.0*thx0).or. & (thy1.eq.-1.0*thy0))) bounceflg2 = .true. if ((thx1.eq.-1.0*thx0).or.(thy1.eq.-1.0*thy0)) bounceflg = .true. c ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************** * Neutrons passing through SEVENTH Al window * 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 alwindow(en0, x0, y0, z0, thx0, thy0, & talx1,talx2,taly1,taly2,talz7, en1, x1, y1, z1, & thx1,thy1,alflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (alflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(alflg1)).or. & ((en0.ne.en1).and.(.not.(alflg1))) ) & write(6,*) ' Edif/alflg Al 7',alflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif Al 7 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.1300) write(6,*) 'Al 7 ',pltot if (dabs(thtot).gt.1.0e9) then write(6,*) ' Al out thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ctemp c if (x1.lt.0) write(6,*) ' Al7 EastDn z1 ',z1,', pathlen ',pltot c if (x1.gt.0) write(6,*) ' Al7 WestUp z1 ',z1,', pathlen ',pltot ctemp ******************************************************************* * Air gap (10cm) before Detector * * 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,tairz6, 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 if (airflg1) beamscflg = .true. if ( ((en0.eq.en1).and.(airflg1)).or. & ((en0.ne.en1).and.(.not.(airflg1))) ) & write(6,*) ' Edif/airflg Al 3',airflg1,en0,en1 c if (en0.ne.en1) write(6,*) ' edif AIR 6 ',airflg1,alflg1, c & beamscflg,tar1flg,tar2flg pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.1300) write(6,*) 'gap 6 ',pltot if (pltot.gt.1300) goto 100 if (x1**2+y1**2.gt.detrad**2) goto 100 Ndetin = Ndetin + wgt DETAVGflux = DETAVGflux + wgt/pi/detrad**2/nevts if (x1**2+y1**2.le.0.625**2) & DETPKflux = DETPKflux + wgt/pi/0.625**2/nevts 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 if (((dabs(x1).lt.1d0).and.(dabs(y1).lt.1d0)).and. & ((dabs(thx1).lt.1d-1).and.(dabs(thy1).lt.1d-1))) & zdet = z1 capfldet = capfldet + sqrt(25.0)/rsqrt(en1) ******************************************************************** ctemp c if (x1.lt.0) write(6,2341) z1,pltot,thtot c if (x1.gt.0) write(6,2342) z1,pltot,thtot ctemp ******************************************************************* c Apply the polarization product -KG lambda=(2.*pi*h_par)/rsqrt(2*m_n*en1) c if((lambda.lt.2.).or.(lambda.gt.14.)) then c write(6,*)'lambda = ', lambda, ' out of range (2,14)Ang' c endif c PA=0.27873+3.5461*(lambda/10.)-6.2591*(lambda/10.)**2 c &+4.5036*(lambda/10.)**3-1.1912*(lambda/10.)**4 c Central PA. See Elog Sim. post 30 cc----------------------------------------------- c Divide into 3 regions: west, middle and east if (PAflg) then if(x1.lt.(-0.95)) then !East PA=-2.8274+47.729*(lambda/10.)-243.93*(lambda/10.)**2 & +614.37*(lambda/10.)**3-845.19*(lambda/10.)**4 & +649.46*(lambda/10.)**5-262.23*(lambda/10.)**6 & +43.397*(lambda/10.)**7 endif if((x1.ge.(-0.95)).and.(x1.le.0.95)) then !Middle PA=0.27873+3.5461*(lambda/10.)-6.2591*(lambda/10.)**2 & +4.5036*(lambda/10.)**3-1.1912*(lambda/10.)**4 endif if(x1.gt.0.95) then !West PA=0.27757+3.3155*(lambda/10.)-5.3255*(lambda/10.)**2 & +3.6243*(lambda/10.)**3-0.94529*(lambda/10.)**4 endif else PA=1 endif cc----------------------------------------------- c if((PA.lt.0.).or.(PA.gt.1.02)) then c write(6,*)'PA = ',PA,' out of range (0,1)' c endif c c c if(vect(6).gt.PA) goto 100 c call hf2(9999,real(lambda),real(PA),wgt) ******************************************************************** *************** End of Beamline ********************************** * c write(6,*) gy*mag_B_beam*480.6*0.1/rsqrt(2.*en1/m_n) , c & gy*mag_B_beam*pltot*0.1/rsqrt(2.*en1/m_n) ! mrad c & , pltot,thtot if((xstart*x1).lt.0.) Nxo=Nxo+wgt !count cross-over neutrons -KG c if((xstart*x1).lt.0.) Nxo=Nxo+1 !-rcm c count bounced neutrons if (bounceflg) numbounce = numbounce + wgt if (bounceflg2) numbounce2 = numbounce2 + wgt if (bounceflg3) numbounce3 = numbounce3 + wgt ******************************************************************** call hf1(3001,real(x1),wgt) call hf1(3002,real(y1),wgt) call hf2(4444,real(x1),real(y1),wgt) call hf1(3003,real(thx1),wgt) call hf1(3004,real(thy1),wgt) call hf1(3005,real(en1),wgt) c lambda=(2.*pi*h_par)/rsqrt(2*m_n*en1) !done in PA application -KG call hf1(4005,real(lambda),wgt) c call hf1(3006,real(pltot),wgt) c call hf1(3007,real(thtot),wgt) ************Rotations and Histograms for ALL and Non-Scattered Neutrons**************** if(x1.ge.0.) then ! x>0 = West (=old right) call hf1(3010,real(pltot),wgt) call hf1(3011,real(thtot),wgt) call hf1(3013,real(en1),wgt) call hf2(6020,real(enstart),real(thtot),wgt) call hf2(6021,real(pltot),real(thtot),wgt) if ((xtest.lt.0).and.(x1.gt.0)) then call hf1(6035, real(pltot),wgt) call hf1(6036, real(thtot),wgt) call hf2(6040,real(enstart),real(thtot),wgt) call hf2(6041,real(pltot),real(thtot),wgt) end if if (.not.((tar1flg.or.tar2flg).and.(vect(5).gt.wgt2))) then ! alternate weighting rall = rall + PA*dsin(thtot*1d-3) !rcm PAsin rallsq = rallsq + (PA*dsin(thtot*1d-3))**2 nrall = nrall + 1 nrallPA = nrallPA + PA nrallPasq = nrallPasq + PA**2 rallavg = rall/nrall rallPAavg = nrallPA/nrall end if rsrot = rsrot + wgt*PA*(dsin(thtot*1d-3)) !rcm PAsin rsrotsq = rsrotsq + wgt*(PA*dsin(thtot*1d-3))**2 nr = nr +wgt nrsq = nrsq + wgt**2 nrPA = nrPA + wgt*PA nrPAsq = nrPAsq + wgt*PA**2 rPAavg = nrPA/nr rsrotavg = rsrot/nr !rcm if (.not.((beamscflg.or.tar1flg).or.tar2flg)) then ! nonscattered rsrotnsc = rsrotnsc + wgt*PA*dsin(thtot*1d-3) rsrotnscsq = rsrotnscsq + wgt*(PA*dsin(thtot*1d-3))**2 nrnsc = nrnsc + wgt nrnscsq = nrnscsq + wgt**2 nrnscPA = nrnscPA + wgt*PA nrnscPAsq = nrnscPAsq + wgt*PA**2 rsnscPAavg = nrnscPA/nrnsc rsrotnscavg = rsrotnsc/nrnsc call hf1(6013, real(pltot),wgt) call hf1(6014, real(thtot),wgt) call hf1(3017,real(en1),wgt) c call hf2(6022,real(entar2),real(thtot),wgt) c call hf2(6023,real(pltot),real(thtot),wgt) c call hf2(6026,real(pltot),real(en1),wgt) rplns = rplns + pltot rplnssq = rplnssq + pltot**2 if (en1.ne.enstart) write(6,*) 'W nonscat delE ',en1-enstart, & beamscflg,tar1flg,tar2flg endif else ! x<0 = East (=old left) call hf1(3008,real(pltot),wgt) call hf1(3009,real(thtot),wgt) call hf1(3012,real(en1),wgt) call hf2(6022,real(enstart),real(thtot),1.0) ! not trajectory dependent -- no wgt call hf2(6023,real(pltot),real(thtot),wgt) if ((xtest.gt.0).and.(x1.lt.0)) then call hf1(6037, real(pltot),wgt) call hf1(6038, real(thtot),wgt) call hf2(6042,real(enstart),real(thtot),1.0) ! not trajectory dependent -- no wgt call hf2(6043,real(pltot),real(thtot),wgt) end if if (.not.((tar1flg.or.tar2flg).and.(vect(5).gt.wgt2))) then ! alternate weighting lall = lall + PA*dsin(thtot*1d-3) lallsq = lallsq + (PA*dsin(thtot*1d-3))**2 nlall = nlall + 1 nlallPA = nlallPA + PA nlallPAsq =nlallPAsq + PA**2 lallavg = lall/nlall lallPAavg = nlallPA/nlall end if lsrot = lsrot + wgt*PA*dsin(thtot*1d-3) !rcm PAsin lsrotsq = lsrotsq + wgt*(PA*dsin(thtot*1d-3))**2 nl = nl + wgt nlsq = nlsq + wgt**2 nlPA = nlPA + wgt*PA nlPAsq = nlPAsq + wgt*PA**2 lPAavg = nlPA/nl lsrotavg = lsrot/nl if (.not.((beamscflg.or.tar1flg).or.tar2flg)) then !rcm PAsin nonscattered lsrotnsc = lsrotnsc + wgt*PA*dsin(thtot*1d-3) lsrotnscsq = lsrotnscsq + wgt*(PA*dsin(thtot*1d-3))**2 nlnsc = nlnsc + wgt nlnscsq = nlnscsq + wgt**2 nlnscPA = nlnscPA + wgt*PA nlnscPAsq = nlnscPAsq + wgt*PA**2 lsnscPAavg = nlnscPA/nlnsc lsrotnscavg = lsrotnsc/nlnsc c call hf1(6009, real(pltot),wgt) c call hf1(6010, real(thtot),wgt) call hf1(6011, real(pltot),wgt) call hf1(6012, real(thtot),wgt) call hf1(3018,real(en1),wgt) c call hf2(6024,real(entar2),real(thtot),wgt) c call hf2(6025,real(pltot),real(thtot),wgt) c call hf2(6027,real(pltot),real(en1),wgt) lplns = lplns + wgt*pltot lplnssq = lplnssq + wgt*pltot**2 if (en1.ne.enstart) write(6,*) 'E nonscat delE ',en1-enstart, & beamscflg,tar1flg,tar2flg endif endif call hf1(7005, real(wgt),1.0) if (wgt.lt.0) write(6,*) ' end loop wgt<0 ',wgt, sig, en1 if (enflg1) numendif1 = numendif1 + wgt if (tar1flg) numtar1flg = numtar1flg + wgt if (enflg2) numendif2 = numendif2 + wgt if (tar2flg) numtar2flg = numtar2flg + wgt if (en1.ne.enstart) numendifdet = numendifdet +wgt if (beamscflg) numbeamsc = numbeamsc + wgt ****************Rotations and Histograms for Scattered neutrons**************** if ((beamscflg.or.tar1flg).or.tar2flg) then ! scattered NEdet=NEdet+wgt call hf1(7006, real(wgt),1.0) c call hf1(6003, real(pltot),wgt) c call hf1(6004, real(thtot),wgt) if (x1.lt.0.) then ! x<0 = East (=old left) call hf2(6031,real(enstart),real(thtot),wgt) call hf1(6007, real(pltot),wgt) call hf1(6008, real(thtot),wgt) if ((xtest.gt.0).and.(x1.lt.0)) then call hf1(6017, real(pltot),wgt) call hf1(6018, real(thtot),wgt) end if call hf1(3015,real(en1),wgt) c call hf2(6024,real(enstart),real(thtot),wgt) c call hf2(6025,real(pltot),real(thtot),wgt) c call hf2(6027,real(pltot),real(en1),wgt) if (.not.((tar1flg.or.tar2flg).and.(vect(5).gt.wgt2))) then ! alternate weighting lsc = lsc + PA*dsin(thtot*1d-3) lscsq = lscsq + (PA*dsin(thtot*1d-3))**2 nlsc = nlsc + 1 nlscPA = nlscPA + PA nlscPAsq = nlscPAsq + PA**2 lscPAavg = nlscPA/nlsc lscavg = lsc/nlsc end if if (tar2flg) then ! scat in East Dn and on E side still call hf1(7007, real(wgt),1.0) lsrottsc = lsrottsc + wgt*PA*dsin(thtot*1d-3) !rcm PAsin lsrottscsq = lsrottscsq + wgt*(PA*dsin(thtot*1d-3))**2 nltsraw = nltsraw + 1 nlts = nlts + wgt nltssq = nltssq + wgt**2 nltsPA = nltsPA + wgt*PA nltsPAsq = nltsPAsq + wgt*PA**2 ltsPAavg = nltsPA/nlts lsrottscavg = lsrottsc/nlts lwde = lwde + wgt*(enstart-en1) lwdesq = lwdesq + wgt*(enstart-en1)**2 lpl = lpl + wgt*pltot lplsq = lplsq + wgt*pltot**2 elseif (tar1flg) then ! scat in West Up i.e. a XTS (cross over tar scat) lsrotxsc = lsrotxsc + wgt*PA*dsin(thtot*1d-3) !rcm PAsin lsrotxscsq = lsrotxscsq + wgt*(PA*dsin(thtot*1d-3))**2 nlxsraw = nlxsraw + 1 nlxs = nlxs + wgt nlxssq = nlxssq + wgt**2 nlxsPA = nlxsPA + wgt*PA nlxsPAsq = nlxsPAsq + wgt*PA**2 lxsPAavg = nlxsPA/nlxs lsrotxscavg = lsrotxsc/nlxs ctemp if (wgt.eq.1.0) write(6,*) ' tar2 wgt=1 ' else ! no target scat, just Al or Gas scat if (.not.(beamscflg)) write(6,*) ' prob with beamscflg ls' lsrotgsc = lsrotgsc + wgt*PA*dsin(thtot*1d-3) !rcm PAsin lsrotgscsq = lsrotgscsq + wgt*(PA*dsin(thtot*1d-3))**2 nlgsraw = nlgsraw + 1 nlgs = nlgs + wgt nlgssq = nlgssq + wgt**2 nlgsPA = nlgsPA + wgt*PA nlgsPAsq = nlgsPAsq + wgt*PA**2 lgsPAavg = nlgsPA/nlgs lsrotgscavg = lsrotgsc/nlgs ctemp if (wgt.ne.1.0) write(6,*) ' tar2 wgt',wgt,tar1flg,tar2flg,xtar1, c & xtar2,x1 end if c call hf1(3020,real(enstart-en1),wgt) call hf1(3021,real(enstart-en1),wgt) c call hf2(8011,real(en1),real(thx1),wgt) c call hf2(8012,real(enstart-en1),real(thx1),wgt) lsde = lsde + wgt*(enstart-en1) lsdesq = lsdesq + wgt*(enstart-en1)**2 ndel = ndel + wgt ndelsq = ndelsq + wgt**2 c write(6,*) ' 1W eloss tar1 tar2 ', entar1-entar2 else ! x>0 = West (=old right) call hf2(6030,real(enstart),real(thtot),wgt) call hf1(6005, real(pltot),wgt) call hf1(6006, real(thtot),wgt) if ((xtest.lt.0).and.(x1.gt.0)) then call hf1(6015, real(pltot),wgt) call hf1(6016, real(thtot),wgt) end if call hf1(3016,real(en1),wgt) c call hf2(6022,real(enstart),real(thtot),wgt) c call hf2(6023,real(pltot),real(thtot),wgt) c call hf2(6026,real(pltot),real(en1),wgt) if (.not.((tar1flg.or.tar2flg).and.(vect(5).gt.wgt2))) then ! alternate weighting rsc = rsc + PA*dsin(thtot*1d-3) rscsq = rscsq + (PA*dsin(thtot*1d-3))**2 nrsc = nrsc + 1 nrscPA = nrscPA + PA nrscPAsq = nrscPAsq + PA**2 rscPAavg = nrscPA/nrsc rscavg = rsc/nrsc end if if (tar1flg) then ! scat West Up and still in West call hf1(7007, real(wgt),1.0) !rcm PAsin rsrottsc = rsrottsc + wgt*PA*dsin(thtot*1d-3) rsrottscsq = rsrottscsq + wgt*(PA*dsin(thtot*1d-3))**2 nrtsraw = nrtsraw + 1 nrts = nrts + wgt nrtssq = nrtssq + wgt**2 nrtsPA = nrtsPA + wgt*PA nrtsPAsq = nrtsPAsq + wgt*PA**2 rtsPAavg = nrtsPA/nrts rsrottscavg = rsrottsc/nrts rwde = rwde + wgt*(enstart-en1) rwdesq = rwdesq + wgt*(enstart-en1)**2 rpl = rpl + wgt*pltot rplsq = rplsq + wgt*pltot**2 else if (tar2flg) then ! scat East Dn so an XTS rsrotxsc = rsrotxsc + wgt*PA*dsin(thtot*1d-3) rsrotxscsq = rsrotxscsq + wgt*(PA*dsin(thtot*1d-3))**2 nrxsraw = nrxsraw + 1 nrxs = nrxs + wgt nrxssq = nrxssq + wgt**2 nrxsPA = nrxsPA + wgt*PA nrxsPAsq = nrxsPAsq + wgt*PA**2 rxsPAavg = nrxsPA/nrxs rsrotxscavg = rsrotxsc/nrxs ctemp if (wgt.eq.1.0) write(6,*) ' tar1 wgt=1 ' else if (.not.(beamscflg)) write(6,*) ' prob with beamscflg rs' rsrotgsc = rsrotgsc + wgt*PA*dsin(thtot*1d-3) !rcm PAsin rsrotgscsq = rsrotgscsq + wgt*(PA*dsin(thtot*1d-3))**2 nrgsraw = nrgsraw + 1 nrgs = nrgs + wgt nrgssq = nrgssq + wgt**2 nrgsPA = nrgsPA + wgt*PA nrgsPAsq = nrgsPAsq + wgt*PA**2 rgsPAavg = nrgsPA/nrgs rsrotgscavg = rsrotgsc/nrgs ctemp if (wgt.ne.1.0) write(6,*) ' tar1 wgt',wgt,tar1flg,tar2flg,xtar1, c & xtar2,x1 end if c call hf1(3020,real(enstart-en1),wgt) call hf1(3022,real(enstart-en1),wgt) c call hf2(8011,real(en1),real(thx1),wgt) c call hf2(8012,real(enstart-en1),real(thx1),wgt) rsde = rsde + wgt*(enstart-en1) rsdesq = rsdesq + wgt*(enstart-en1)**2 nder = nder + wgt ndersq = ndersq + wgt**2 c write(6,*) ' 2E eloss tar1 tar2 ', entar1-entar2 endif endif ******************************************************************** * alambda = 2*pi*h_par/rsqrt(2*m_n*en1) c ispec = int(alambda/0.016) ispec = int(alambda/0.030) if (ispec.gt.1000) then write(6,*) ' ispec:',ispec,' set to 1000, alambda:', & alambda ispec = 1000 end if spec(ispec,1) = alambda spec(ispec,2) = spec(ispec,2) + 1 ******************************************************************** if ((xtar1*xtar2.gt.0.0).and. ! xtar1 = x-pos before UPtar & ((x1*xtar1.lt.0.0).or.(x1*xtar2.lt.0.0)) ) then c write(6,*) xtar1, xtar2, x1 dncross = dncross + wgt end if c if (abs(thtot).gt.300.0) then c write(6,941) xtest,x1,enstart,en1,pltot,thtot c endif 941 format(1x,1pe9.2,2x,1pe9.2,2x,1pe9.2,2x,1pe9.2,2x, & 1pe12.5,2x,1pe9.2) if ((xtest.lt.0).and.(x1.gt.0)) numE2W = numE2W + wgt if ((xtest.gt.0).and.(x1.lt.0)) numW2E = numW2E + wgt 34 format(i12,2x,a6) ctemp c if (x1.lt.0) write(6,*) 'End EastDn ',z1,pltot ,thtot c if (x1.gt.0) write(6,*) 'End WestUp ',z1,pltot,thtot ctemp 100 do j=1,10 lim = float(j)*float(nevts)/10. if(float(i).eq.lim) then write(*,34) i, 'events' c write(*,*) 'rsrot=',rsrot,' rsrotsq=',rsrotsq,' nrPA=',nrPA !-KG !rcm c write(*,*) 'lsrot=',lsrot,' lsrotsq=',lsrotsq,' nlPa=',nlPA !-KG !rcm endif enddo 101 continue do j=1, 1000 cbec write(71,*) j, spec(j,1), spec(j,2) enddo c write(*,*) 'calling nhe3' call nhe3(spec) print * print *, " through loop, writing to ", filename write(6,*) ' num of tar1flg ',numtar1flg write(6,*) ' num of endif 1',numendif1 write(3,*) ' num of tar1flg ',numtar1flg write(3,*) ' num of endif 1',numendif1 write(6,*) ' num of tar2flg ',numtar2flg write(6,*) ' num of endif 2',numendif2 write(3,*) ' num of tar2flg ',numtar2flg write(3,*) ' num of endif 2',numendif2 write(6,*) ' num of beamscflg ',numbeamsc write(3,*) ' num of beamscflg ',numbeamsc write(6,*) ' num Tar scat West ',nrts write(6,*) ' num Tar scat East ',nlts write(3,*) ' num Tar scat West ',nrts write(3,*) ' num Tar scat East ',nlts write(6,*) ' num Gas scat West ',nrgs write(6,*) ' num Gas scat East ',nlgs write(3,*) ' num Gas scat West ',nrgs write(3,*) ' num Gas scat East ',nlgs write(6,*) ' num of endet.ne.enstart ',numendifdet write(3,*) ' num of endet.ne.enstart ',numendifdet write(6,*) ' num E to W at det ',numE2W write(6,*) ' num W to E at det ',numW2E write(3,*) ' num E to W at det ',numE2W write(3,*) ' num W to E at det ',numW2E print * if (hegastar) then write(6,*) ' 4He gas target ', ndenhegas write(3,*) ' 4He gas target ', ndenhegas elseif (airtar) then write(6,*) ' Air target ', ndenair write(3,*) ' Air target ', ndenair else write(6,*) ' Liquid He target; Temp= ',T,' K ' write(3,*) ' Liquid He target; Temp= ',T,' K ' endif print * write(*,*) 'input ',nevts write(*,*) 'PSM output ',NPSMout,NPSMout/NPSMout write(*,*) 'Air gap 1 out ',Nairout1,Nairout1/NPSMout write(*,*) 'Al window 1 out ',Nalout1,Nalout1/NPSMout write(*,*) 'Inguide Output ',Nguideout,Nguideout/NPSMout ! rcm write(*,*) 'Al window 2 out ',Nalout2,Nalout2/NPSMout ! rcm write(*,*) 'Input Coil output ',Nincoilout,Nincoilout/NPSMout write(*,*) 'WE target input ',NtarWEin,NtarWEin/NPSMout write(*,*) 'Up West target IN ',Ntar1,Ntar1/NPSMout write(*,*) 'Up West target OUT ',Ntar1out,Ntar1out/NPSMout write(*,*) 'Up East empty IN ',Nemp1in,Nemp1in/NPSMout write(*,*) 'Up East empty Out ',Nemp1out,Nemp1out/NPSMout write(*,*) 'Dn West empty IN ',Nemp2in,Nemp2in/NPSMout write(*,*) 'DN West empty Out ',Nemp2out,Nemp2out/NPSMout write(*,*) 'Dn East target IN ',Ntar2,Ntar2/NPSMout write(*,*) 'Dn East target OUT ',Ntar2out,Ntar2out/NPSMout write(*,*) 'E-W Target Outputs ',NTargetout,NTargetout/NPSMout write(*,*) 'Output Coil output ',NOutcoilout,NOutcoilout/NPSMout write(*,*) 'ASM output (geom.) ', dble(cot), cot/NPSMout write(*,*) 'ASM out scattered ',NEdet,NEdet/NPSMout write(*,*) 'Detector input ',Ndetin,Ndetin/NPSMout write(*,*) write(*,5551) 0.25/xmax/ymax,INPKflux,0.25/xmax/ymax/INPKflux write(*,5552) PSMAVGflux,PSMPKflux,PSMAVGflux/PSMPKflux write(*,5553) INCOAVGflux,INCOPKflux,INCOAVGflux/INCOPKflux write(*,5554) TARAVGflux,TARPKflux,TARAVGFlux/TARPKFlux write(*,5555) OUTCOAVGflux,OUTCOPKflux,OUTCOAVGflux/OUTCOPKflux write(*,5556) DETAVGflux,DETPKflux,DETAVGFlux/DETPKFlux write(*,*) write(*,5560) INCOAVGFlux/PSMAVGflux,INCOPKFlux/PSMPKflux, & INCOPKFlux/PSMPKflux/capflpsm*nPSMout*capflinc/nincoilout write(*,5561) OUTCOAVGFlux/PSMAVGflux,OUTCOPKFlux/PSMPKflux, & OUTCOPKFlux/PSMPKflux/capflpsm*nPSMout*capflout/NOutcoilout write(*,5562) DETAVGFlux/PSMAVGflux,DETPKFlux/PSMPKflux, & DETPKFlux/PSMPKflux/capflpsm*nPSMout*capflout/NOutcoilout write(*,*) write(3,5560) INCOAVGFlux/PSMAVGflux,INCOPKFlux/PSMPKflux, & INCOPKFlux/PSMPKflux/capflpsm*nPSMout*capflinc/nincoilout write(3,5561) OUTCOAVGFlux/PSMAVGflux,OUTCOPKFlux/PSMPKflux, & OUTCOPKFlux/PSMPKflux/capflpsm*nPSMout*capflout/NOutcoilout write(3,5562) DETAVGFlux/PSMAVGflux,DETPKFlux/PSMPKflux, & DETPKFlux/PSMPKflux/capflpsm*nPSMout*capflout/NOutcoilout write(3,*) 5560 format(1x,' Incoil/PSM av',1x,1pe10.3,', pk',1x,1pe10.3, & ', CAPpk',1x,1pe10.3) 5561 format(1x,' Outcoil/PSM av',1x,1pe10.3,', pk',1x,1pe10.3, & ', CAPpk',1x,1pe10.3) 5562 format(1x,' Detector/PSM av',1x,1pe10.3,', pk',1x,1pe10.3, & ', CAPpk',1x,1pe10.3) write(*,*) ' Cap Flux ratio PSM ', capflpsm/nPSMout write(*,*) ' Cap Flux ratio Incoil ', capflinc/nincoilout write(*,*) ' Cap Flux ratio Outcoil ', capflout/NOutcoilout write(*,*) ' Cap Flux ratio Detector ', capfldet/ndetin write(*,*) write(*,5651) xmax,ymax write(*,5652) xin,yin if (sptmflg) write(*,5659) ix1,ix2,ix3,ix4,iy1,iy2 if (.not.(sptmflg)) write(*,5653) ix1,ix4,iy1,iy2 write(*,5654) tarx1,tarx2,tarx3,tarx4 write(*,5655) tary1,tary2 write(*,5656) ox1,ox2,ox3,ox4 write(*,5657) oy1,oy2 write(*,5658) detrad write(*,*) 5551 format(1x,' Input Flux av',1x,1pe10.3,', pk',1x,1pe10.3, & ', av/pk',1x,1pe10.3) 5552 format(1x,' PSM Flux av',1x,1pe10.3,', pk',1x,1pe10.3, & ', av/pk',1x,1pe10.3) 5553 format(1x,' InCoil Flux av',1x,1pe10.3,', pk',1x,1pe10.3, & ', av/pk',1x,1pe10.3) 5554 format(1x,' Target Flux av',1x,1pe10.3,', pk',1x,1pe10.3, & ', av/pk',1x,1pe10.3) 5555 format(1x,' Outcoil Flux av',1x,1pe10.3,', pk',1x,1pe10.3, & ', av/pk',1x,1pe10.3) 5556 format(1x,' Detector Flux av',1x,1pe10.3,', pk',1x,1pe10.3, & ', av/pk',1x,1pe10.3) 5651 format(1x,' Input Dim (cm)',1x,1f5.2,1x,1f5.2) 5652 format(1x,' PSM Dim (cm)',1x,1f5.2,1x,1f5.2) 5653 format(1x,' Incoil Dim (cm)',1x,1f5.2,1x,1f5.2,1x,1f5.2, & 1x,1f5.2) 5659 format(1x,' Incoil Dim (cm)',1x,1f5.2,1x,1f5.2,1x,1f5.2, & 1x,1f5.2,1x,1f5.2,1x,1f5.2) 5654 format(1x,' Target Dim X (cm)',1x,1f5.2,1x,1f5.2,1x,1f5.2, & 1x,1f5.2) 5655 format(1x,' Target Dim Y (cm)',1x,1f5.2,1x,1f5.2) 5656 format(1x,' Outcol Dim X (cm)',1x,1f5.2,1x,1f5.2,1x,1f5.2, & 1x,1f5.2) 5657 format(1x,' Outcoil Dim Y (cm)',1x,1f5.2,1x,1f5.2) 5658 format(1x,' Detector radius (cm)',1x,1f5.2) write(3,*) 'input ',nevts write(3,*) 'PSM output ',NPSMout,NPSMout/NPSMout write(3,*) 'Air gap 1 out ',Nairout1,Nairout1/NPSMout write(3,*) 'Al window 1 out ',Nalout1,Nalout1/NPSMout write(3,*) 'Inguide Output ',Nguideout,Nguideout/NPSMout ! rcm write(3,*) 'Al window 2 out ',Nalout2,Nalout2/NPSMout ! rcm write(3,*) 'Input Coil output ',Nincoilout,Nincoilout/NPSMout write(3,*) 'WE target input ',NtarWEin,NtarWEin/NPSMout write(3,*) 'Up West target IN ',Ntar1,Ntar1/NPSMout write(3,*) 'Up West target OUT ',Ntar1out,Ntar1out/NPSMout write(3,*) 'Up East empty IN ',Nemp1in,Nemp1in/NPSMout write(3,*) 'Up East empty Out ',Nemp1out,Nemp1out/NPSMout write(3,*) 'Dn West empty IN ',Nemp2in,Nemp2in/NPSMout write(3,*) 'DN West empty Out ',Nemp2out,Nemp2out/NPSMout write(3,*) 'Dn East target IN ',Ntar2,Ntar2/NPSMout write(3,*) 'Dn East target OUT ',Ntar2out,Ntar2out/NPSMout write(3,*) 'E-W Target Outputs ',NTargetout,NTargetout/NPSMout write(3,*) 'Output Coil output ',NOutcoilout,NOutcoilout/NPSMout write(3,*) 'ASM output (geom.) ', dble(cot), cot/NPSMout write(3,*) 'ASM out scattered ',NEdet,NEdet/NPSMout write(3,*) 'Detector input ',Ndetin,Ndetin/NPSMout write(3,*) write(3,5551) 0.25/xmax/ymax,INPKflux,0.25/xmax/ymax/INPKflux write(3,5552) PSMAVGflux,PSMPKflux,PSMAVGflux/PSMPKflux write(3,5553) INCOAVGflux,INCOPKflux,INCOAVGflux/INCOPKflux write(3,5554) TARAVGflux,TARPKflux,TARAVGFlux/TARPKFlux write(3,5555) OUTCOAVGflux,OUTCOPKflux,OUTCOAVGflux/OUTCOPKflux write(3,5556) DETAVGflux,DETPKflux,DETAVGFlux/DETPKFlux write(3,*) write(3,*) ' Incoil/PSM av',INCOAVGFlux/PSMAVGflux, & ', pk',INCOPKFlux/PSMPKflux write(3,*) ' Outcoil/PSM av',OUTCOAVGFlux/PSMAVGflux, & ', pk',OUTCOPKFlux/PSMPKflux write(3,*) ' Detector/PSM av',DETAVGFlux/PSMAVGflux, & ', pk',DETPKFlux/PSMPKflux write(3,*) write(3,*) ' Cap Flux ratio PSM ', capflpsm/nPSMout write(3,*) ' Cap Flux ratio Incoil ', capflinc/nincoilout write(3,*) ' Cap Flux ratio Outcoil ', capflout/NOutcoilout write(3,*) ' Cap Flux ratio Detector ', capfldet/ndetin write(3,*) write(3,5651) xmax,ymax write(3,5652) xin,yin if (sptmflg) write(3,5659) ix1,ix2,ix3,ix4,iy1,iy2 if (.not.(sptmflg)) write(3,5653) ix1,ix4,iy1,iy2 write(3,5654) tarx1,tarx2,tarx3,tarx4 write(3,5655) tary1,tary2 write(3,5656) ox1,ox2,ox3,ox4 write(3,5657) oy1,oy2 write(3,5658) detrad write(3,*) write(*,*) 'Number scattered in Al window 1 ',Nalscat1 write(3,*) 'Number scattered in Al window 1 ',Nalscat1 write(*,*) 'Number scattered in Al window 2 ',Nalscat2 write(3,*) 'Number scattered in Al window 2 ',Nalscat2 !rcm write(*,*) 'Number scattered in Air gap 1 ',Nairscat1 write(3,*) 'Number scattered in Air gap 1 ',Nairscat1 write(*,*) 'Number scattered in UpTar ',Nenscat1 write(3,*) 'Number scattered in UpTar ',Nenscat1 write(*,*) 'Number scattered in DnTar ',Nenscat2 write(3,*) 'Number scattered in DnTar ',Nenscat2 write(*,*) 'Number scattered in EmptyTar Up ',Nempty1 write(*,*) 'Number scattered in EmptyTar Dn ',Nempty2 write(3,*) 'Number scattered in EmptyTar Up ',Nempty1 write(3,*) 'Number scattered in EmptyTar Dn ',Nempty2 c write(*,*) 'Number Theta scattered in Target 1 ',Nthscat1 c write(3,*) 'Number Theta scattered in Target 1 ',Nthscat1 c write(*,*) 'Number Theta scattered in Target 2 ',Nthscat2 c write(3,*) 'Number Theta scattered in Target 2 ',Nthscat2 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 Adding some printouts at end of run -- bec print * Call Date_and_time(date,time,zone,datetime) write(6,*) ' End time: ', date, time write(6,*) ' HBOOK filename: ',outfname write(6,*) c*********averages -- Note: Now with weights, nr, nl, etc are all sums of c******* weights, see logbook (bec) fsw = nrts/nr ufsw = fsw*rsqrt(1.0/dble(nrtsraw)+1.0/dble(nr)) fse = nlts/nl ufse = fse*rsqrt(1.0/dble(nltsraw)+1.0/dble(nl)) fxw = nrxs/nr ufxw = fxw*rsqrt(1.0/dble(nrxsraw)+1.0/dble(nr)) fxe = nlxs/nl ufxe = fxe*rsqrt(1.0/dble(nlxsraw)+1.0/dble(nl)) write(6,755) Nenscat1/Ntar1 write(6,756) Nenscat2/Ntar2 write(6,757) Ntar12/Ntar1 write(6,758) cot/float(Nevts) write(6,759) fsw,ufsw write(6,760) fse,ufse write(6,761) nrgs/nr,nrgs/nr*dsqrt(1/nrgsraw+1/nr) write(6,762) nlgs/nl,nlgs/nl*dsqrt(1/nlgsraw+1/nl) write(6,765) fxw,ufxw write(6,766) fxe,ufxe write(6,767) nrnsc/nr,nrnsc/nr*dsqrt(1/nrnsc+1/nr) write(6,768) nlnsc/nl,nlnsc/nl*dsqrt(1/nlnsc+1/nl) write(6,*) write(6,*) ' StDev All West ', & rsqrt((rsrotsq-rsrot**2/nr)/nr)*1.0d3 write(6,*) ' St Err All West ', & rsqrt((rsrotsq-rsrot**2/nr)/nr)*rSqrt(nrsq)/nr*1.0d3 !rcm c write(6,*) ' nrPA ',nrPA , ' nrPAsq ',nrPAsq c write(6,*) 'rsqrt(nrPAsq)/nrPA ', rSqrt(nrPAsq)/nrPA c write(6,*) ' StDev NoScat West ', & rsqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc)*1.0d3 write(6,*) ' St Err NoScat West ', & rsqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc)* & rSqrt(nrnscsq)/nrnsc*1.0d3 c write(3,*) write(3,*) ' StDev XOvTar Scat West ', & rsqrt((rsrotxscsq-rsrotxsc**2/nrxs)/nrxs)*1.0d3 write(3,*) ' St Err XOvTar Scat West ', & rsqrt((rsrotxscsq-rsrotxsc**2/nrxs)/nrxs)*1.0d3 & *rSqrt(nrxssq)/nrxs c write(3,*) ' nrxs ',nrxs , ' nrxssq ',nrxssq c write(3,*) 'rsqrt(nrxssq)/nrxs ', rSqrt(nrxssq)/nrxs if ((rsrottscsq-rsrottsc**2/nrts)/nrts.gt.0) then c write(6,*) write(6,*) ' StDev Tar Scat West ', & rsqrt((rsrottscsq-rsrottsc**2/nrts)/nrts)*1.0d3 write(6,*) ' St Err Tar West ', & rsqrt((rsrottscsq-rsrottsc**2/nrts)/nrts)*1.0d3 & *rSqrt(nrtssq)/nrts c write(6,*) ' nrts ',nrts , ' nrtssq ',nrtssq c write(6,*) 'rsqrt(nrtssq)/nrts ', rSqrt(nrtssq)/nrts c write(6,*) end if if ((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs.gt.0) then write(6,*) ' StDev Gas Scat West ', & rsqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs)*1.0d3 write(6,*) ' St Err Gas West ', & rsqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs)*1.0d3 & *rSqrt(nrgssq)/nrgs c write(6,*) ' nrgs ',nrgs , ' nrgssq ',nrgssq c write(6,*) 'rsqrt(nrgssq)/nrgs ', rSqrt(nrgssq)/nrgs end if write(6,*) write(6,*) write(6,713) nr, nl write(6,714) nr/nl write(6,717) nr-nl,rsqrt(nlsq+nrsq) write(6,*) write(6,735) NguideoutW,NguideoutE write(6,736) NincoiloutW,NincoiloutE write(6,737) NoutcoiloutW,NoutcoiloutE write(6,*) write(6,770) nrnsc , nlnsc write(6,*) write(6,700) zdet, mag_B,mag_B_beam, bzgradnum write(6,*) write(6,*) ' ********************** Average Rotation Values &***********************' write(6,*) argrsrot = (rsrotsq-rsrot**2/nr)/nr argrsrotPA = (nrPAsq-nrPA**2/nr)/nr c if (argrsrot.lt.0) then c write(6,*) 'argrsrot',argrsrot c argrsrot =0 c endif c if (argrsrotPA.lt.0) then c write(6,*) 'argrsrotPA',argrsrotPA c argrsrotPA =0 c endif uncrsrotavg = rsqrt(argrsrot)*rsqrt(nrsq)/nr uncrPAavg = rsqrt(argrsrotPA)*rsqrt(nrsq)/nr arglsrot = (lsrotsq-lsrot**2/nl)/nl arglsrotPA = (nlPAsq-nlPA**2/nl)/nl c if (arglsrot.lt.0) then c write(6,*) 'arglsrot',arglsrot c arglsrot =0 c endif c if (arglsrotPA.lt.0) then c write(6,*) 'arglsrotPA',arglsrotPA c arglsrotPA =0 c endif unclsrotavg = rsqrt(arglsrot)*rsqrt(nlsq)/nl unclPAavg = rsqrt(arglsrotPA)*rsqrt(nlsq)/nl write(6,704) real(rPAavg*100),dasin(rsrotavg/rPAavg)*1d3, !rcm & rsqrt((uncrsrotavg/rPAavg)**2+ & (uncrPAavg*rsrotavg/rPAavg**2)**2)/ & rsqrt(1-(rsrotavg/rPAavg)**2)*1d3 write(6,703) real(lPAavg*100),dasin(lsrotavg/lPAavg)*1d3, & rsqrt((unclsrotavg/lPAavg)**2+ & (unclPAavg*lsrotavg/lPAavg**2)**2)/ & rsqrt(1-(lsrotavg/lPAavg)**2)*1d3 write(6,711) dasin(rsrotavg/rPAavg)*1d3- & dasin(lsrotavg/lPAavg)*1d3, & rsqrt( ((uncrsrotavg/rPAavg)**2+ & (uncrPAavg*rsrotavg/rPAavg**2)**2)/ & (1-(rsrotavg/rPAavg)**2) + ((unclsrotavg/lPAavg)**2+ & (unclPAavg*lsrotavg/lPAavg**2)**2)/ & (1-(lsrotavg/lPAavg)**2))*1d3 c write(6,704) rsrot/nr, c & rsqrt((rsrotsq-rsrot**2/nr)/nr)*rSqrt(nrsq)/nr c write(6,703) lsrot/nl, c & rsqrt((lsrotsq-lsrot**2/nl)/nl)*rSqrt(nlsq)/nl c write(6,711) rsrot/nr-lsrot/nl, c & rsqrt( (lsrotsq-lsrot**2/nl)/nl*nlsq/nl**2 c & + (rsrotsq-rsrot**2/nr)/nr*nrsq/nr**2 ) c if ((nrall.gt.1).and.(nlall.gt.1)) then c write(6,*) c write(6,804) rall/float(nrall), ! alternate weighting method c & rsqrt((rallsq-rall**2/nrall)/nrall/(nrall-1)) c write(6,803) lall/nlall, c & rsqrt((lallsq-lall**2/nlall)/nlall/(nlall-1)) c write(6,811) rall/nrall-lall/nlall, c & rsqrt( (lallsq-lall**2/nlall)/nlall/(nlall-1) c & + (rallsq-rall**2/nrall)/nrall/(nlall-1) ) c else c write(6,*) 'No alternate weighting data ' c end if uncrallavg = rsqrt((rallsq-rall**2/nrall)/ !rcm PAsin & (nrall-1)/nrall) uncrallPAavg = rsqrt((nrallPAsq - nrallPA**2/nrall)/ & (nrall-1)/nrall) unclallavg = rsqrt((lallsq-lall**2/nlall)/ & (nlall-1)/nlall) unclallPAavg = rsqrt((nlallPAsq - nlallPA**2/nlall)/ & (nlall-1)/nlall) write(6,*) if ((nrall.gt.1).and.(nlall.gt.1)) then write(6,804) real(rallPAavg*100),dasin(rallavg/rallPAavg)*1d3, ! alternate weighting method & rsqrt( (uncrallavg/rallPAavg)**2+ & (uncrallPAavg*rallavg/rallPAavg**2)**2 )/ & rsqrt(1-(rallavg/rallPAavg)**2)*1d3 write(6,803) real(lallPAavg*100),dasin(lallavg/lallPAavg)*1d3, & rsqrt( (unclallavg/lallPAavg)**2+ & (unclallPAavg*lallavg/lallPAavg**2)**2 )/ & rsqrt(1-(lallavg/lallPAavg)**2)*1d3 write(6,811) dasin(rallavg/rallPAavg)*1d3- & dasin(lallavg/lallPAavg)*1d3, & rsqrt( ((uncrallavg/rallPAavg)**2+ & (uncrallPAavg*rallavg/rallPAavg**2)**2)/ & (1-(rallavg/rallPAavg)**2) + ((unclallavg/lallPAavg)**2+ & (unclallPAavg*lallavg/lallPAavg**2)**2)/ & (1-(lallavg/lallPAavg)**2))*1d3 else write(6,*) 'No alternate weighting data ' end if write(6,*) argrsnsc = (rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc argrsnscPA = (nrnscPAsq-nrnscPA**2/nrnsc)/nrnsc c if (argrsnsc.lt.0) then c write(6,*) 'argrsnsc',argrsnsc c argrsnsc =0 c endif c if (argrsnscPA.lt.0) then c write(6,*) 'argrsnscPA',argrsnscPA c argrsnscPA =0 c endif uncrsnsc = rsqrt(argrsnsc)*rsqrt(nrnscsq)/nrnsc uncrsnscPA = rsqrt(argrsnscPA)*rsqrt(nrnscsq)/nrnsc arglsnsc = (lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc arglsnscPA = (nlnscPAsq-nlnscPA**2/nlnsc)/nlnsc c if (arglsnsc.lt.0) then c write(6,*) 'arglsnsc',arglsnsc c arglsnsc =0 c endif c if (arglsnscPA.lt.0) then c write(6,*) 'arglsnscPA',arglsnscPA c arglsnscPA =0 c endif unclsnsc = rsqrt(arglsnsc)*rsqrt(nlnscsq)/nlnsc unclsnscPA = rsqrt(arglsnsc)*rsqrt(nlnscsq)/nlnsc write(6,726) real(rsnscPAavg*100), & dasin(rsrotnscavg/rsnscPAavg)*1d3, & rsqrt((uncrsnsc/rsnscPAavg)**2+ & (uncrsnscPA*rsrotnscavg/rsnscPAavg**2)**2)/ & rsqrt(1-(rsrotnscavg/rsnscPAavg)**2)*1d3 write(6,725) real(lsnscPAavg*100), & dasin(lsrotnscavg/lsnscPAavg)*1d3, & rsqrt((unclsnsc/lsnscPAavg)**2+ !rcm & (unclsnscPA*lsrotavg/lsnscPAavg**2)**2)/ & rsqrt(1-(lsrotnscavg/lsnscPAavg)**2)*1d3 write(6,732)dasin(rsrotnscavg/rsnscPAavg)*1d3- & dasin(lsrotnscavg/lsnscPAavg)*1d3, & rsqrt( ((uncrsnsc/rsnscPAavg)**2+ & (uncrsnscPA*rsrotnscavg/rsnscPAavg**2)**2)/ & (1-(rsrotnscavg/rsnscPAavg)**2) + ((unclsnsc/lsnscPAavg)**2+ & (unclsnscPA*lsrotnscavg/lsnscPAavg**2)**2)/ & (1-(lsrotnscavg/lsnscPAavg)**2) )*1d3 c write(6,726) rsrotnsc/nrnsc, !rcm c & rsqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc) c & *rSqrt(nrnscsq)/nrnsc c write(6,725) lsrotnsc/nlnsc, c & rsqrt((lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc) c & *rSqrt(nlnscsq)/nlnsc c write(6,732) rsrotnsc/nrnsc - lsrotnsc/nlnsc, c & rsqrt( (lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc c & *nlnscsq/nlnsc**2 c & + (rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc c & *nrnscsq/nrnsc**2) write(6,*) if ((nrsc.gt.1).and.(nlsc.gt.1)) then !rcm uncrscavg = rsqrt((rscsq-rsc**2/nrsc)/ & (nrsc-1)/nrsc) uncrscPAavg = rsqrt((nrscPAsq - nrscPA**2/nrsc)/ & (nrsc-1)/nrsc) unclscavg = rsqrt((lscsq-lsc**2/nlsc)/ & (nlsc-1)/nlsc) unclscPAavg = rsqrt((nlscPAsq - nlscPA**2/nlsc)/ & (nlsc-1)/nlsc) write(6,806) real(rscPAavg*100),dasin(rscavg/rscPAavg)*1d3, ! alternate weighting method & rsqrt( (uncrscavg/rscPAavg)**2+ & (uncrscPAavg*rscavg/rscPAavg**2)**2 )/ & rsqrt(1-(rscavg/rscPAavg)**2)*1d3 write(6,805) real(lscPAavg*100),dasin(lscavg/lscPAavg)*1d3, & rsqrt( (unclscavg/lscPAavg)**2+ & (unclscPAavg*lscavg/lscPAavg**2)**2 )/ & rsqrt(1-(lscavg/lscPAavg)**2)*1d3 write(6,812) dasin(rscavg/rscPAavg)*1d3-dasin(lscavg/lscPAavg)*1d3, & rsqrt( ((uncrscavg/rscPAavg)**2+(uncrscPAavg*rscavg/ & rscPAavg**2)**2)/(1-(rscavg/rscPAavg)**2) + & ((unclscavg/lscPAavg)**2+(unclscPAavg*lscavg/ & lscPAavg**2)**2)/(1-(lscavg/lscPAavg)**2) )*1d3 else write(6,*) ' No alternate weighting scatters' end if write(6,*) if ((nlts.gt.0).and.(nrts.gt.0)) then argrtsc = (rsrottscsq-rsrottsc**2/nrts)/nrts argrtscPA = (nrtsPAsq-nrtsPA**2/nrts)/nrts uncrstsc = rsqrt(argrtsc)*rsqrt(nrtssq)/nrts uncrstscPA = rsqrt(argrtscPA)*rsqrt(nrtssq)/nrts argltsc = (lsrottscsq-lsrottsc**2/nlts)/nlts argltscPA = (nltsPAsq-nltsPA**2/nlts)/nlts unclstsc = rsqrt(argltsc)*rsqrt(nltssq)/nlts unclstscPA = rsqrt(argltscPA)*rsqrt(nltssq)/nlts write(6,706) real(rtsPAavg*100),dasin(rsrottscavg/rtsPAavg)*1d3, & rsqrt((uncrstsc/rtsPAavg)**2+ & (uncrstscPA*rsrottscavg/rtsPAavg**2)**2)/ & rsqrt(1-(rsrottscavg/rtsPAavg)**2)*1d3 write(6,705) real(ltsPAavg*100),dasin(lsrottscavg/ltsPAavg)*1d3, & rsqrt((unclstsc/ltsPAavg)**2+ & (unclstscPA*lsrottscavg/ltsPAavg**2)**2)/ & rsqrt(1-(lsrottscavg/ltsPAavg)**2)*1d3 write(6,712) dasin(rsrottscavg/rtsPAavg)*1d3- & dasin(lsrottscavg/ltsPAavg)*1d3, & rsqrt( ((uncrstsc/rtsPAavg)**2+(uncrstscPA*rsrottscavg/ & rtsPAavg**2)**2)/(1-(rsrottscavg/rtsPAavg)**2) + & ((unclstsc/ltsPAavg)**2+(unclstscPA*lsrottscavg/ & ltsPAavg**2)**2)/(1-(lsrottscavg/ltsPAavg)**2) )*1d3 end if c write(6,706) rsrottsc/nrts, c & rsqrt((rsrottscsq-rsrottsc**2/nrts)/nrts) c & *rSqrt(nrtssq)/nrts c write(6,705) lsrottsc/nlts, c & rsqrt((lsrottscsq-lsrottsc**2/nlts)/nlts) c & *rSqrt(nltssq)/nrts c write(6,712) rsrottsc/nrts - lsrottsc/nlts , c & rsqrt( (lsrottscsq-lsrottsc**2/nlts)/nlts c & *nltssq/nlts**2 c & + (rsrottscsq-rsrottsc**2/nrts)/nrts c & *nrtssq/nrts**2 ) write(6,*) if ((nlgs.gt.0).and.(nrgs.gt.0)) then argrgsc = (rsrotgscsq-rsrotgsc**2/nrgs)/nrgs argrgscPA = (nrgsPAsq-nrgsPA**2/nrgs)/nrgs uncrsgsc = rsqrt(argrgsc)* rsqrt(nrgssq)/nrgs uncrsgscPA = rsqrt(argrgscPA)*rsqrt(nrgssq)/nrgs arglgsc = (lsrotgscsq-lsrotgsc**2/nlgs)/nlgs arglgscPA = (nlgsPAsq-nlgsPA**2/nlgs)/nlgs unclsgsc = rsqrt(arglgsc)* rsqrt(nlgssq)/nlgs unclsgscPA = rsqrt(arglgscPA)*rsqrt(nlgssq)/nlgs write(6,7006) real(rgsPAavg*100),dasin(rsrotgscavg/rgsPAavg)*1d3, & rsqrt((uncrsgsc/rgsPAavg)**2+ & (uncrsgscPA*rsrotgscavg/rgsPAavg**2)**2)/ & rsqrt(1-(rsrotgscavg/rgsPAavg)**2)*1d3 write(6,7005) real(lgsPAavg*100),dasin(lsrotgscavg/lgsPAavg)*1d3, & rsqrt((unclsgsc/lgsPAavg)**2+ & (unclsgscPA*lsrotgscavg/lgsPAavg**2)**2)/ & rsqrt(1-(lsrotgscavg/lgsPAavg)**2)*1d3 write(6,7012) dasin(rsrotgscavg/rgsPAavg)*1d3- & dasin(lsrotgscavg/lgsPAavg)*1d3, & rsqrt( ((uncrsgsc/rgsPAavg)**2+(uncrsgscPA*rsrotgscavg/ & rgsPAavg**2)**2)/(1-(rsrotgscavg/rgsPAavg)**2) + & ((unclsgsc/lgsPAavg)**2+(unclsgscPA*lsrotgscavg/ & lgsPAavg**2)**2)/(1-(lsrotgscavg/lgsPAavg)**2) )*1d3 end if write(6,*) if ((nlxs.gt.0).and.(nrxs.gt.0)) then argrxsc = (rsrotxscsq-rsrotxsc**2/nrxs)/nrxs argrxscPA = (nrxsPAsq-nrxsPA**2/nrxs)/nrxs uncrsxsc = rsqrt(argrxsc)* rsqrt(nrxssq)/nrxs uncrsxscPA = rsqrt(argrxscPA)*rsqrt(nrxssq)/nrxs arglxsc = (lsrotxscsq-lsrotxsc**2/nlxs)/nlxs arglxscPA = (nlxsPAsq-nlxsPA**2/nlxs)/nlxs unclsxsc = rsqrt(arglxsc)* rsqrt(nlxssq)/nlxs unclsxscPA = rsqrt(arglxscPA)*rsqrt(nlxssq)/nlxs write(6,7008) real(rxsPAavg*100),dasin(rsrotxscavg/rxsPAavg)*1d3, & rsqrt((uncrsxsc/rxsPAavg)**2+ & (uncrsxscPA*rsrotxscavg/rxsPAavg**2)**2)/ & rsqrt(1-(rsrotxscavg/rxsPAavg)**2)*1d3 write(6,7007) real(lxsPAavg*100),dasin(lsrotxscavg/lxsPAavg)*1d3, & rsqrt((unclsxsc/lxsPAavg)**2+ & (unclsxscPA*lsrotxscavg/lxsPAavg**2)**2)/ & rsqrt(1-(lsrotxscavg/lxsPAavg)**2)*1d3 write(6,7013) dasin(rsrotxscavg/rxsPAavg)*1d3- & dasin(lsrotxscavg/lxsPAavg)*1d3, & rsqrt( ((uncrsxsc/rxsPAavg)**2+(uncrsxscPA*rsrotxscavg/ & rxsPAavg**2)**2)/(1-(rsrotxscavg/rxsPAavg)**2) + & ((unclsxsc/lxsPAavg)**2+(unclsxscPA*lsrotxscavg/ & lxsPAavg**2)**2)/(1-(lsrotxscavg/lxsPAavg)**2) )*1d3 end if c write(6,7006) rsrotgsc/nrgs, c & rsqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs) c & *rSqrt(nrgssq)/nrgs c write(6,7005) lsrotgsc/nlgs c & rsqrt((lsrotgscsq-lsrotgsc**2/nlgs)/nlgs) c & *rSqrt(nlgssq)/nrgs c write(6,7012) rsrotgsc/nrgs - lsrotgsc/nlgs , c & rsqrt( (lsrotgscsq-lsrotgsc**2/nlgs)/nlgs c & *nlgssq/nlgs**2 c & + (rsrotgscsq-rsrotgsc**2/nrgs)/nrgs c & *nrgssq/nrgs**2 ) cccccccc estimated PNC angle scaling TScat and Xscat into average of ALL if (nlxs.gt.0) then rotnsw = dasin(rsrotnscavg/rsnscPAavg)*1d3 urotnsw = rsqrt((uncrsnsc/rsnscPAavg)**2+ & (uncrsnscPA*rsrotnscavg/rsnscPAavg**2)**2)/ & rsqrt(1-(rsrotnscavg/rsnscPAavg)**2)*1d3 rotnse =dasin(lsrotnscavg/lsnscPAavg)*1d3 urotnse = rsqrt((unclsnsc/lsnscPAavg)**2+ !rcm & (unclsnscPA*lsrotnscavg/lsnscPAavg**2)**2)/ & rsqrt(1-(lsrotnscavg/lsnscPAavg)**2)*1d3 rotsw = dasin(rsrottscavg/rtsPAavg)*1d3 urotsw = rsqrt((uncrstsc/rtsPAavg)**2+ & (uncrstscPA*rsrottscavg/rtsPAavg**2)**2)/ & rsqrt(1-(rsrottscavg/rtsPAavg)**2)*1d3 rotse = dasin(lsrottscavg/ltsPAavg)*1d3 urotse = rsqrt((unclstsc/ltsPAavg)**2+ & (unclstscPA*lsrottscavg/ltsPAavg**2)**2)/ & rsqrt(1-(lsrottscavg/ltsPAavg)**2)*1d3 rotxw = dasin(rsrotxscavg/rxsPAavg)*1d3 urotxw = rsqrt((uncrsxsc/rxsPAavg)**2+ & (uncrsxscPA*rsrotxscavg/rxsPAavg**2)**2)/ & rsqrt(1-(rsrotxscavg/rxsPAavg)**2)*1d3 rotxe = dasin(lsrotxscavg/lxsPAavg)*1d3 urotxe = rsqrt((unclsxsc/lxsPAavg)**2+ & (unclsxscPA*lsrotxscavg/lxsPAavg**2)**2)/ & rsqrt(1-(lsrotxscavg/lxsPAavg)**2)*1d3 write(6,*) write(6,*) write(6,889) 0.5*( (fse-fsw)*0.5*(rotnsw+rotnse) + & fsw*rotsw-fse*rotse ), & 0.5*rsqrt( ((0.5*(rotnsw+rotnse)-rotse)*ufse)**2 + & (-(0.5*(rotnsw+rotnse)+rotsw)*ufsw)**2 + & (0.5*(fse-fsw)*urotnse)**2 + (0.5*(fse-fsw)*urotnsw)**2 + & (fsw*urotsw)**2 + (fse*urotse)**2 ) write(6,*) write(6,888) 0.5*( (fse-fsw+fxe-fxw)*0.5*(rotnsw+rotnse) + & fsw*rotsw-fse*rotse + fxw*rotxw-fxe*rotxe ), & 0.5*rsqrt( ((0.5*(rotnsw+rotnse)-rotse)*ufse)**2 + & (-(0.5*(rotnsw+rotnse)+rotsw)*ufsw)**2 + & (0.5*(fse-fsw+fxe-fxw)*urotnse)**2 + & (0.5*(fse-fsw+fxe-fxw)*urotnsw)**2 + & (fsw*urotsw)**2 +(fse*urotse)**2 + & (fxw*urotxw)**2 +(fxe*urotxe)**2 + & (rotxw*ufxw)**2 +(rotxe*ufxe)**2 ) write(6,*) write(6,890) (fse-fsw+fxe-fxw)*0.5*(rotnsw+rotnse), & rsqrt( (ufsw**2+ufse**2+ufxw**2+ufxe**2)* & (0.5*(rotnsw+rotnse))**2 & + ((fse-fsw+fxe-fxw)*0.5)**2*(urotnsw**2+urotnse**2) ) write(6,893) (fse-fsw)*0.5*(rotnsw+rotnse), & rsqrt( (ufsw**2+ufse**2)*(0.5*(rotnsw+rotnse))**2 & + ((fse-fsw)*0.5)**2*(urotnsw**2+urotnse**2)**2 ) write(6,891) fsw*rotsw-fse*rotse , & rsqrt( (ufsw*rotsw)**2 + (fsw*urotsw)**2 & + (ufse*rotse)**2 + (fse*urotse)**2 ) write(6,892) fxw*rotxw-fxe*rotxe , & rsqrt( (ufxw*rotxw)**2 + (fxw*urotxw)**2 & + (ufxe*rotxe)**2 + (fxe*urotxe)**2 ) write(6,*) end if write(6,*) ' ***************************************************', &'******************' write(6,*) if (nder.gt.0) then write(6,702) rsde/nder, & rsqrt((rsdesq-rsde**2/nder)/nder)*rSqrt(ndersq)/nder write(6,701) lsde/ndel, & rsqrt((lsdesq-lsde**2/ndel)/ndel)*rSqrt(ndelsq)/ndel end if if ((nrts.gt.1).and.(nlts.gt.1)) then write(6,802) rwde/nrts, & rsqrt((rwdesq-rwde**2/nrts)/nrts)*rsqrt(nrtssq)/nrts write(6,801) lwde/nlts, & rsqrt((lwdesq-lwde**2/nlts)/nlts)*rsqrt(nltssq)/nlts write(6,822) rpl/nrts, & rsqrt((rplsq-rpl**2/nrts)/nrts)*rsqrt(nrtssq)/nrts write(6,821) lpl/nlts, & rsqrt((lplsq-lpl**2/nlts)/nlts)*rsqrt(nltssq)/nlts end if write(6,*) write(6,*) ' NonScat avg pathlen (cm) West ',rplns/nrnsc,' +- ', & rsqrt((rplnssq-rplns**2/nrnsc)/nrnsc)*rSqrt(nrnscsq)/nrnsc write(6,*) ' NonScat avg pathlen (cm) East ',lplns/nlnsc,' +- ', & rsqrt((lplnssq-lplns**2/nlnsc)/nlnsc)*rSqrt(nlnscsq)/nlnsc write(6,*) write(6,715) nrts, nlts write(6,716) nrts/nlts write(6,718) nrts-nlts,rSqrt(nltssq+nrtssq) write(6,782) nrtsraw,nltsraw write(6,*) write(6,7015) nrgs, nlgs write(6,7016) nrgs/nlgs write(6,7018) nrgs-nlgs,rSqrt(nlgssq+nrgssq) write(6,7082) nrgsraw,nlgsraw write(6,*) write(6,7020) nrxs, nlxs write(6,7021) nrxs/nlxs write(6,7022) nrxs-nlxs,rsqrt(nlxssq+nrxssq) write(6,7023) nrxsraw,nlxsraw write(3,*) if ((nrsc.gt.0).and.(nlsc.gt.0)) then write(6,815) nrsc, nlsc write(6,816) nrsc/nlsc write(6,818) nrsc-nlsc,sqrt(float(nlsc+nrsc)) else write(6,*) ' NO scattering so limited output written ' endif c 700 format(3x,' Rotations for Incoil Ouput to Outcoil input 114.6~cm with c & B(mGauss)=',1f9.4) 700 format(3x,' PSMout to Det ',f9.4,' ~cm ', & ' Bt(mG)=',1f6.2,'Bb(mG)=',1f6.2,' BzGradNum ',i2) 701 format(3x,' Avg E loss Det East Dn All Scat (meV)',1pe11.4,' +- ' & ,1pe11.4) 702 format(3x,' Avg E loss Det West Up All Scat (meV)',1pe11.4,' +- ' & ,1pe11.4) 801 format(3x,' Avg E loss Det EDn Target Scat (meV)',1pe11.4,' +- ' & ,1pe11.4) 802 format(3x,' Avg E loss Det WUp Target Scat (meV)',1pe11.4,' +- ' & ,1pe11.4) 821 format(3x,' Avg PathLength Det EDn Tar Scat (meV)',1pe15.8,' +- ' & ,1pe11.4) 822 format(3x,' Avg PathLength Det WUp Tar Scat (meV)',1pe15.8,' +- ' & ,1pe11.4) 703 format(3x,' PA ',f5.1,' East Dn (mrad) ',1pe13.6, & ' +- ',1pe9.2) 704 format(3x,' PA ',f5.1,' West Up (mrad) ',1pe13.6, & ' +- ',1pe9.2) 803 format(3x,' PA ',f5.1,' Rot East Dn (alt wgt) (mrad) ',1pe13.6, & ' +- ',1pe9.2) 804 format(3x,' PA ',f5.1,' Rot West Up (alt wgt) (mrad) ',1pe13.6, & ' +- ',1pe9.2) 705 format(3x,' PA ',f5.1,' Target Scat East Dn (mrad) ',1pe13.6, & ' +- ',1pe9.2) 706 format(3x,' PA ',f5.1,' Target Scat West Up (mrad) ',1pe13.6, & ' +- ',1pe9.2) 805 format(3x,' PA ',f5.1,' Scat East Dn (alt wgt) (mrad) ',1pe13.6, & ' +- ',1pe9.2) 806 format(3x,' PA ',f5.1,' Scat West Up (alt wgt) (mrad) ',1pe13.6, & ' +- ',1pe9.2) 7005 format(3x,' PA ',f5.1,' Gas Scat East Dn (mrad) ',1pe13.6, & ' +- ',1pe9.2) 7006 format(3x,' PA ',f5.1,' Gas Scat West Up (mrad) ',1pe13.6, & ' +- ',1pe9.2) 7007 format(3x,' PA ',f5.1,' Xover Tar Scat East Dn (mrad) ',1pe13.6, & ' +- ',1pe9.2) 7008 format(3x,' PA ',f5.1,' Xover Tar Scat West Up (mrad) ',1pe13.6, & ' +- ',1pe9.2) 725 format(3x,' PA ',f5.1,' NOScat East Dn (mrad) ',1pe13.6, & ' +- ',1pe9.2) 726 format(3x,' PA ',f5.1,' NOScat West Up (mrad) ',1pe13.6, & ' +- ',1pe9.2) c 707 format(3x,' Avg. Rot Scat East Dn E<0.35meV (mrad) ',1pe11.4, c & ' +- ',1pe11.4) c 708 format(3x,' Avg. Rot Scat West Up E<0.35meV (mrad) ',1pe11.4, c & ' +- ',1pe11.4) c 709 format(3x,' Frac of rotation from scat East Dn ',1pe11.4, c & ', oflow En ',1pe11.4) c 710 format(3x,' Frac of rotation from scat West Up ',1pe11.4, c & ', oflow En ',1pe11.4) 711 format(3x,' UpW-DnE (mrad) ',1pe13.6, & ' +- ',1pe9.2) 811 format(3x,' UpW-DnE (alt wgt) (mrad) ',1pe13.6, & ' +- ',1pe9.2) 712 format(3x,' UpW-DnE Tar scat (mrad) ',1pe13.6, & ' +- ',1pe9.2) 812 format(3x,' UpW-DnE Scat (alt wgt) (mrad) ',1pe13.6, & ' +- ',1pe9.2) 7012 format(3x,' UpW-DnE Gas scat (mrad) ',1pe13.6, & ' +- ',1pe9.2) 7013 format(3x,' UpW-DnE Xov Tar Scat (mrad) ',1pe13.6, & ' +- ',1pe9.2) 888 format(3x,' Estimated PNC angle from TS and XTS (mrad) ',1pe13.6, & ' +- ',1pe9.2) 889 format(3x,' Estimated PNC angle from TS only (mrad) ',1pe13.6, & ' +- ',1pe9.2) 890 format(3x,' (fse-fsw+fxe-fxw)*Theta-ns-avg (mrad) ',1pe13.6, & ' +- ',1pe9.2) 893 format(3x,' (fse-fsw)*Theta-ns-avg (mrad) ',1pe13.6, & ' +- ',1pe9.2) 891 format(3x,' fsw*thsw-fse*thse (mrad) ',1pe13.6, & ' +- ',1pe9.2) 892 format(3x,' fxw*thxw-fxe*thxe (mrad) ',1pe13.6, & ' +- ',1pe9.2) 732 format(3x,' UpW-DnE NOScat (mrad) ',1pe13.6, & ' +- ',1pe9.2) 713 format(3x,' # in DetW Up',1pe13.6,2x,', # in DetE Dn ',1pe13.6) 714 format(3x,' Ratio # DetWUp/DetEDn ,'1pe13.6) 717 format(3x,' Det Counts UpW-DnE ',1pe11.4,' +- ',1pe11.4) 715 format(3x,' # Tarscat in DetW Up ',1pe11.4,2x, & ', # Tarscat in Det E Dn ',1pe11.4) 716 format(3x,' Ratio # Tarscat DetWUp/DetEDn ,'1pe11.4) 718 format(3x,' TarScat Det Counts UpW-DnE ',1pe11.4,' +- ',1pe11.4) 815 format(3x,' AltWgt # Tarscat in DetW Up ',1pe11.4,2x, & ', # Tarscat in Det E Dn ',1pe11.4) 816 format(3x,' AltWgt Ratio # Tarscat DetWUp/DetEDn ,'1pe11.4) 818 format(3x,' AltWgt TarScat Det Counts UpW-DnE ',1pe11.4,' +- ' & ,1pe11.4) 7015 format(2x,'# Gas/Al scat DetW Up ',1pe11.4,2x, & ',# Gas/Al scat DetE Dn ',1pe11.4) 7016 format(3x,' Ratio # Gas/Al Scat DetWUp/DetEDn ,'1pe11.4) 7018 format(3x,' G/Al Scat Det Counts UpW-DnE ',1pe11.4,' +- ', & 1pe11.4) 7020 format(2x,'# XovTarScat scat DetW Up ',1pe11.4,2x, & ',# XovTarScat scat DetE Dn ',1pe11.4) 7021 format(3x,' Ratio # XovTarScat DetWUp/DetEDn ,'1pe11.4) 7022 format(3x,' XovTarScat Det Counts UpW-DnE ',1pe11.4,' +- ', & 1pe11.4) 719 format(3x,'# scat w/ E<0.35meV at det West: ',1pe11.4, & ' East: ',1pe11.4) 755 format(1x,' Frac # of scatters with Ef.ne.Ei tar1 ',1pe11.4) 756 format(1x,' Frac # of scatters with Ef.ne.Ei tar2 ',1pe11.4) 757 format(1x,' Frac num of tar1 entering tar2 ',1pe11.4) 758 format(1x,' Frac # of total reaching det ',1pe11.4) 759 format(1x,' Frac det hits TarScat W ',1pe11.4,' +- ',1pe11.4) 760 format(1x,' Frac det hits TarScat E ',1pe11.4,' +- ',1pe11.4) 761 format(1x,' Frac det hits G/Al Scat W ',1pe11.4,' +- ',1pe11.4) 762 format(1x,' Frac det hits G/Al Scat E ',1pe11.4,' +- ',1pe11.4) 763 format(1x,' Frac det hits TXall W ',1pe11.4,' +- ',1pe11.4) 764 format(1x,' Frac det hits TXall E ',1pe11.4,' +- ',1pe11.4) 765 format(1x,' Frac det hits TXScat W ',1pe11.4,' +- ',1pe11.4) 766 format(1x,' Frac det hits TXScat E ',1pe11.4,' +- ',1pe11.4) 767 format(1x,' Frac det hits NO Scat W ',1pe11.4,' +- ',1pe11.4) 768 format(1x,' Frac det hits NO Scat E ',1pe11.4,' +- ',1pe11.4) 770 format(3x,' num NO scat -- det W: ',1pe11.4, & ' det E: ',1pe11.4) 780 format(3x,'Total No. of cross-over neutrons: ',1pe11.4) 782 format(3x,'Unwt # Tar scat DetW Up ',1pe11.4,2x, & ', Unwt # Tar scat Det E Dn ',1pe11.4) 7082 format(3x,'UnWt # Gas scat DetW Up ',1pe11.4,2x, & ', Unwt # Gas scat Det E Dn ',1pe11.4) 7023 format(3x,'UnWt # XovTarScat DetW Up ',1pe11.4,2x, & ', Unwt # XovTarScat Det E Dn ',1pe11.4) 735 format(3x,' Guide output W ',1pe11.4,' Guide output E ',1pe11.4) 736 format(3x,' Incoil output W ',1pe11.4,' Incoil output E ',1pe11.4) 737 format(3x,' Outcoil output W ',1pe11.4, & ' Outcoil output E ',1pe11.4) c write(3,*) ' End time: ', date, time write(3,*) ' HBOOK filename: ',outfname write(3,*) c*********averages -- Note: Now with weights, nr, nl, etc are all sums of c******* weights, see logbook (bec) write(3,755) Nenscat1/Ntar1 write(3,756) Nenscat2/Ntar2 write(3,757) Ntar12/Ntar1 write(3,758) cot/float(Nevts) write(3,759) fsw,ufsw write(3,760) fse,ufse write(3,761) nrgs/nr,nrgs/nr*dsqrt(1/nrgsraw+1/nr) write(3,762) nlgs/nl,nlgs/nl*dsqrt(1/nlgsraw+1/nl) write(3,765) fxw,ufxw write(3,766) fxe,ufxe write(3,767) nrnsc/nr,nrnsc/nr*dsqrt(1/nrnsc+1/nr) write(3,768) nlnsc/nl,nlnsc/nl*dsqrt(1/nlnsc+1/nl) write(3,*) write(3,*) ' StDev All West ', & rsqrt((rsrotsq-rsrot**2/nr)/nr)*1d3 write(3,*) ' St Err All West ', & rsqrt((rsrotsq-rsrot**2/nr)/nr)*rSqrt(nrsq)/nr*1d3 c write(3,*) ' nr ',nr , ' nrsq ',nrsq c write(3,*) 'rsqrt(nrsq)/nr ', rSqrt(nrsq)/nr c write(3,*) ' StDev NoScat West ', & rsqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc)*1d3 write(3,*) ' St Err NoScat West ', & rsqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc)* & rSqrt(nrnscsq)/nrnsc*1d3 c write(3,*) ' nr ',nr , ' nrsq ',nrsq c write(3,*) 'rsqrt(nrsq)/nr ', rSqrt(nrsq)/nr c if ((rsrottscsq-rsrottsc**2/nrts)/nrts.gt.0) then c write(3,*) write(3,*) ' StDev Tar Scat West ', & rsqrt((rsrottscsq-rsrottsc**2/nrts)/nrts)*1d3 write(3,*) ' St Err Tar Scat West ', & rsqrt((rsrottscsq-rsrottsc**2/nrts)/nrts)*1d3 & *rSqrt(nrtssq)/nrts c write(3,*) ' nrts ',nrts , ' nrtssq ',nrtssq c write(3,*) 'rsqrt(nrtssq)/nrts ', rSqrt(nrtssq)/nrts end if c write(3,*) write(3,*) ' StDev XOvTar Scat West ', & rsqrt((rsrotxscsq-rsrotxsc**2/nrxs)/nrxs)*1d3 write(3,*) ' St Err XOvTar Scat West ', & rsqrt((rsrotxscsq-rsrotxsc**2/nrxs)/nrxs)*1d3 & *rSqrt(nrxssq)/nrxs c write(3,*) ' nrxs ',nrxs , ' nrxssq ',nrxssq c write(3,*) 'rsqrt(nrxssq)/nrxs ', rSqrt(nrxssq)/nrxs if ((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs.gt.0) then c write(3,*) write(3,*) ' StDev Gas Scat West ', & rsqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs)*1d3 write(3,*) ' St Err Gas Scat West ', & rsqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs)*1d3 & *rSqrt(nrgssq)/nrgs c write(3,*) ' nrgs ',nrgs , ' nrgssq ',nrgssq c write(3,*) 'rsqrt(nrgssq)/nrgs ', rSqrt(nrgssq)/nrgs end if write(3,*) write(3,*) write(3,713) nr, nl write(3,714) nr/nl write(3,717) nr-nl,rSqrt(nlsq+nrsq) write(3,*) write(3,735) NguideoutW,NguideoutE write(3,736) NincoiloutW,NincoiloutE write(3,737) NoutcoiloutW,NoutcoiloutE write(3,*) write(3,770) nrnsc , nlnsc write(3,*) write(3,700) zdet, mag_B,mag_B_beam, bzgradnum write(3,*) write(3,*) ' ********************** Average Rotation Values &***********************' write(3,*) argrsrot = (rsrotsq-rsrot**2/nr)/nr argrsrotPA = (nrPAsq-nrPA**2/nr)/nr c if (argrsrot.lt.0) then c write(3,*) 'argrsrot',argrsrot c argrsrot =0 c endif c if (argrsrotPA.lt.0) then c write(3,*) 'argrsrotPA',argrsrotPA c argrsrotPA =0 c endif uncrsrotavg = rsqrt(argrsrot)*rsqrt(nrsq)/nr uncrPAavg = rsqrt(argrsrotPA)*rsqrt(nrsq)/nr arglsrot = (lsrotsq-lsrot**2/nl)/nl arglsrotPA = (nlPAsq-nlPA**2/nl)/nl c if (arglsrot.lt.0) then c write(3,*) 'arglsrot',arglsrot c arglsrot =0 c endif c if (arglsrotPA.lt.0) then c write(3,*) 'arglsrotPA',arglsrotPA c arglsrotPA =0 c endif unclsrotavg = rsqrt(arglsrot)*rsqrt(nlsq)/nl unclPAavg = rsqrt(arglsrotPA)*rsqrt(nlsq)/nl write(3,704) real(rPAavg*100),dasin(rsrotavg/rPAavg)*1d3, !rcm & rsqrt((uncrsrotavg/rPAavg)**2+ & (uncrPAavg*rsrotavg/rPAavg**2)**2)/ & rsqrt(1-(rsrotavg/rPAavg)**2)*1d3 write(3,703) real(lPAavg*100),dasin(lsrotavg/lPAavg)*1d3, & rsqrt((unclsrotavg/lPAavg)**2+ & (unclPAavg*lsrotavg/lPAavg**2)**2)/ & rsqrt(1-(lsrotavg/lPAavg)**2)*1d3 write(3,711) dasin(rsrotavg/rPAavg)*1d3- & dasin(lsrotavg/lPAavg)*1d3, & rsqrt( ((uncrsrotavg/rPAavg)**2+ & (uncrPAavg*rsrotavg/rPAavg**2)**2)/ & (1-(rsrotavg/rPAavg)**2) + ((unclsrotavg/lPAavg)**2+ & (unclPAavg*lsrotavg/lPAavg**2)**2)/ & (1-(lsrotavg/lPAavg)**2))*1d3 c write(3,704) rsrot/nr, c & rsqrt((rsrotsq-rsrot**2/nr)/nr)*rSqrt(nrsq)/nr c write(3,703) lsrot/nl, c & rsqrt((lsrotsq-lsrot**2/nl)/nl)*rSqrt(nlsq)/nl c write(3,711) rsrot/nr-lsrot/nl, c & rsqrt( (lsrotsq-lsrot**2/nl)/nl*nlsq/nl**2 c & + (rsrotsq-rsrot**2/nr)/nr*nrsq/nr**2 ) c if ((nrall.gt.1).and.(nlall.gt.1)) then c write(3,*) c write(3,804) rall/float(nrall), ! alternate weighting method c & rsqrt((rallsq-rall**2/nrall)/nrall/(nrall-1)) c write(3,803) lall/nlall, c & rsqrt((lallsq-lall**2/nlall)/nlall/(nlall-1)) c write(3,811) rall/nrall-lall/nlall, c & rsqrt( (lallsq-lall**2/nlall)/nlall/(nlall-1) c & + (rallsq-rall**2/nrall)/nrall/(nlall-1) ) c else c write(3,*) 'No alternate weighting data ' c end if uncrallavg = rsqrt((rallsq-rall**2/nrall)/ !rcm PAsin & (nrall-1)/nrall) uncrallPAavg = rsqrt((nrallPAsq - nrallPA**2/nrall)/ & (nrall-1)/nrall) unclallavg = rsqrt((lallsq-lall**2/nlall)/ & (nlall-1)/nlall) unclallPAavg = rsqrt((nlallPAsq - nlallPA**2/nlall)/ & (nlall-1)/nlall) write(3,*) if ((nrall.gt.1).and.(nlall.gt.1)) then write(3,804) real(rallPAavg*100),dasin(rallavg/rallPAavg)*1d3, ! alternate weighting method & rsqrt( (uncrallavg/rallPAavg)**2+ & (uncrallPAavg*rallavg/rallPAavg**2)**2 )/ & rsqrt(1-(rallavg/rallPAavg)**2)*1d3 write(3,803) real(lallPAavg*100),dasin(lallavg/lallPAavg)*1d3, & rsqrt( (unclallavg/lallPAavg)**2+ & (unclallPAavg*lallavg/lallPAavg**2)**2 )/ & rsqrt(1-(lallavg/lallPAavg)**2)*1d3 write(3,811) dasin(rallavg/rallPAavg)*1d3- & dasin(lallavg/lallPAavg)*1d3, & rsqrt( ((uncrallavg/rallPAavg)**2+ & (uncrallPAavg*rallavg/rallPAavg**2)**2)/ & (1-(rallavg/rallPAavg)**2) + ((unclallavg/lallPAavg)**2+ & (unclallPAavg*lallavg/lallPAavg**2)**2)/ & (1-(lallavg/lallPAavg)**2))*1d3 else write(3,*) 'No alternate weighting data ' end if write(3,*) argrsnsc = (rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc argrsnscPA = (nrnscPAsq-nrnscPA**2/nrnsc)/nrnsc c if (argrsnsc.lt.0) then c write(3,*) 'argrsnsc',argrsnsc c argrsnsc =0 c endif c if (argrsnscPA.lt.0) then c write(3,*) 'argrsnscPA',argrsnscPA c argrsnscPA =0 c endif uncrsnsc = rsqrt(argrsnsc)*rsqrt(nrnscsq)/nrnsc uncrsnscPA = rsqrt(argrsnscPA)*rsqrt(nrnscsq)/nrnsc arglsnsc = (lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc arglsnscPA = (nlnscPAsq-nlnscPA**2/nlnsc)/nlnsc c if (arglsnsc.lt.0) then c write(3,*) 'arglsnsc',arglsnsc c arglsnsc =0 c endif c if (arglsnscPA.lt.0) then c write(3,*) 'arglsnscPA',arglsnscPA c arglsnscPA =0 c endif unclsnsc = rsqrt(arglsnsc)*rsqrt(nlnscsq)/nlnsc unclsnscPA = rsqrt(arglsnsc)*rsqrt(nlnscsq)/nlnsc write(3,726) real(rsnscPAavg*100), & dasin(rsrotnscavg/rsnscPAavg)*1d3, & rsqrt((uncrsnsc/rsnscPAavg)**2+ & (uncrsnscPA*rsrotnscavg/rsnscPAavg**2)**2)/ & rsqrt(1-(rsrotnscavg/rsnscPAavg)**2)*1d3 write(3,725) real(lsnscPAavg*100), & dasin(lsrotnscavg/lsnscPAavg)*1d3, & rsqrt((unclsnsc/lsnscPAavg)**2+ !rcm & (unclsnscPA*lsrotavg/lsnscPAavg**2)**2)/ & rsqrt(1-(lsrotnscavg/lsnscPAavg)**2)*1d3 write(3,732)dasin(rsrotnscavg/rsnscPAavg)*1d3- & dasin(lsrotnscavg/lsnscPAavg)*1d3, & rsqrt( ((uncrsnsc/rsnscPAavg)**2+ & (uncrsnscPA*rsrotnscavg/rsnscPAavg**2)**2)/ & (1-(rsrotnscavg/rsnscPAavg)**2) + ((unclsnsc/lsnscPAavg)**2+ & (unclsnscPA*lsrotnscavg/lsnscPAavg**2)**2)/ & (1-(lsrotnscavg/lsnscPAavg)**2) )*1d3 c write(3,726) rsrotnsc/nrnsc, !rcm c & rsqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc) c & *rSqrt(nrnscsq)/nrnsc c write(3,725) lsrotnsc/nlnsc, c & rsqrt((lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc) c & *rSqrt(nlnscsq)/nlnsc c write(3,732) rsrotnsc/nrnsc - lsrotnsc/nlnsc, c & rsqrt( (lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc c & *nlnscsq/nlnsc**2 c & + (rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc c & *nrnscsq/nrnsc**2) write(3,*) if ((nrsc.gt.1).and.(nlsc.gt.1)) then !rcm uncrscavg = rsqrt((rscsq-rsc**2/nrsc)/ & (nrsc-1)/nrsc) uncrscPAavg = rsqrt((nrscPAsq - nrscPA**2/nrsc)/ & (nrsc-1)/nrsc) unclscavg = rsqrt((lscsq-lsc**2/nlsc)/ & (nlsc-1)/nlsc) unclscPAavg = rsqrt((nlscPAsq - nlscPA**2/nlsc)/ & (nlsc-1)/nlsc) write(3,806) real(rscPAavg*100),dasin(rscavg/rscPAavg)*1d3, ! alternate weighting method & rsqrt( (uncrscavg/rscPAavg)**2+ & (uncrscPAavg*rscavg/rscPAavg**2)**2 )/ & rsqrt(1-(rscavg/rscPAavg)**2)*1d3 write(3,805) real(lscPAavg*100),dasin(lscavg/lscPAavg)*1d3, & rsqrt( (unclscavg/lscPAavg)**2+ & (unclscPAavg*lscavg/lscPAavg**2)**2 )/ & rsqrt(1-(lscavg/lscPAavg)**2)*1d3 write(3,812) dasin(rscavg/rscPAavg)*1d3-dasin(lscavg/lscPAavg)*1d3, & rsqrt( ((uncrscavg/rscPAavg)**2+(uncrscPAavg*rscavg/ & rscPAavg**2)**2)/(1-(rscavg/rscPAavg)**2) + & ((unclscavg/lscPAavg)**2+(unclscPAavg*lscavg/ & lscPAavg**2)**2)/(1-(lscavg/lscPAavg)**2) )*1d3 else write(3,*) ' No alternate weighting scatters' end if write(3,*) if ((nlts.gt.0).and.(nrts.gt.0)) then argrtsc = (rsrottscsq-rsrottsc**2/nrts)/nrts argrtscPA = (nrtsPAsq-nrtsPA**2/nrts)/nrts uncrstsc = rsqrt(argrtsc)*rsqrt(nrtssq)/nrts uncrstscPA = rsqrt(argrtscPA)*rsqrt(nrtssq)/nrts argltsc = (lsrottscsq-lsrottsc**2/nlts)/nlts argltscPA = (nltsPAsq-nltsPA**2/nlts)/nlts unclstsc = rsqrt(argltsc)*rsqrt(nltssq)/nlts unclstscPA = rsqrt(argltscPA)*rsqrt(nltssq)/nlts write(3,706) real(rtsPAavg*100),dasin(rsrottscavg/rtsPAavg)*1d3, & rsqrt((uncrstsc/rtsPAavg)**2+ & (uncrstscPA*rsrottscavg/rtsPAavg**2)**2)/ & rsqrt(1-(rsrottscavg/rtsPAavg)**2)*1d3 write(3,705) real(ltsPAavg*100),dasin(lsrottscavg/ltsPAavg)*1d3, & rsqrt((unclstsc/ltsPAavg)**2+ & (unclstscPA*lsrottscavg/ltsPAavg**2)**2)/ & rsqrt(1-(lsrottscavg/ltsPAavg)**2)*1d3 write(3,712) dasin(rsrottscavg/rtsPAavg)*1d3- & dasin(lsrottscavg/ltsPAavg)*1d3, & rsqrt( ((uncrstsc/rtsPAavg)**2+(uncrstscPA*rsrottscavg/ & rtsPAavg**2)**2)/(1-(rsrottscavg/rtsPAavg)**2) + & ((unclstsc/ltsPAavg)**2+(unclstscPA*lsrottscavg/ & ltsPAavg**2)**2)/(1-(lsrottscavg/ltsPAavg)**2) )*1d3 end if c write(3,706) rsrottsc/nrts, c & rsqrt((rsrottscsq-rsrottsc**2/nrts)/nrts) c & *rSqrt(nrtssq)/nrts c write(3,705) lsrottsc/nlts, c & rsqrt((lsrottscsq-lsrottsc**2/nlts)/nlts) c & *rSqrt(nltssq)/nrts c write(3,712) rsrottsc/nrts - lsrottsc/nlts , c & rsqrt( (lsrottscsq-lsrottsc**2/nlts)/nlts c & *nltssq/nlts**2 c & + (rsrottscsq-rsrottsc**2/nrts)/nrts c & *nrtssq/nrts**2 ) write(3,*) if ((nlgs.gt.0).and.(nrgs.gt.0)) then argrgsc = (rsrotgscsq-rsrotgsc**2/nrgs)/nrgs argrgscPA = (nrgsPAsq-nrgsPA**2/nrgs)/nrgs uncrsgsc = rsqrt(argrgsc)* rsqrt(nrgssq)/nrgs uncrsgscPA = rsqrt(argrgscPA)*rsqrt(nrgssq)/nrgs arglgsc = (lsrotgscsq-lsrotgsc**2/nlgs)/nlgs arglgscPA = (nlgsPAsq-nlgsPA**2/nlgs)/nlgs unclsgsc = rsqrt(arglgsc)* rsqrt(nlgssq)/nlgs unclsgscPA = rsqrt(arglgscPA)*rsqrt(nlgssq)/nlgs write(3,7006) real(rgsPAavg*100),dasin(rsrotgscavg/rgsPAavg)*1d3, & rsqrt((uncrsgsc/rgsPAavg)**2+ & (uncrsgscPA*rsrotgscavg/rgsPAavg**2)**2)/ & rsqrt(1-(rsrotgscavg/rgsPAavg)**2)*1d3 write(3,7005) real(lgsPAavg*100),dasin(lsrotgscavg/lgsPAavg)*1d3, & rsqrt((unclsgsc/lgsPAavg)**2+ & (unclsgscPA*lsrotgscavg/lgsPAavg**2)**2)/ & rsqrt(1-(lsrotgscavg/lgsPAavg)**2)*1d3 write(3,7012) dasin(rsrotgscavg/rgsPAavg)*1d3- & dasin(lsrotgscavg/lgsPAavg)*1d3, & rsqrt( ((uncrsgsc/rgsPAavg)**2+(uncrsgscPA*rsrotgscavg/ & rgsPAavg**2)**2)/(1-(rsrotgscavg/rgsPAavg)**2) + & ((unclsgsc/lgsPAavg)**2+(unclsgscPA*lsrotgscavg/ & lgsPAavg**2)**2)/(1-(lsrotgscavg/lgsPAavg)**2) )*1d3 end if write(3,*) if ((nlxs.gt.0).and.(nrxs.gt.0)) then argrxsc = (rsrotxscsq-rsrotxsc**2/nrxs)/nrxs argrxscPA = (nrxsPAsq-nrxsPA**2/nrxs)/nrxs uncrsxsc = rsqrt(argrxsc)* rsqrt(nrxssq)/nrxs uncrsxscPA = rsqrt(argrxscPA)*rsqrt(nrxssq)/nrxs arglxsc = (lsrotxscsq-lsrotxsc**2/nlxs)/nlxs arglxscPA = (nlxsPAsq-nlxsPA**2/nlxs)/nlxs unclsxsc = rsqrt(arglxsc)* rsqrt(nlxssq)/nlxs unclsxscPA = rsqrt(arglxscPA)*rsqrt(nlxssq)/nlxs write(3,7008) real(rxsPAavg*100),dasin(rsrotxscavg/rxsPAavg)*1d3, & rsqrt((uncrsxsc/rxsPAavg)**2+ & (uncrsxscPA*rsrotxscavg/rxsPAavg**2)**2)/ & rsqrt(1-(rsrotxscavg/rxsPAavg)**2)*1d3 write(3,7007) real(lxsPAavg*100),dasin(lsrotxscavg/lxsPAavg)*1d3, & rsqrt((unclsxsc/lxsPAavg)**2+ & (unclsxscPA*lsrotxscavg/lxsPAavg**2)**2)/ & rsqrt(1-(lsrotxscavg/lxsPAavg)**2)*1d3 write(3,7013) dasin(rsrotxscavg/rxsPAavg)*1d3- & dasin(lsrotxscavg/lxsPAavg)*1d3, & rsqrt( ((uncrsxsc/rxsPAavg)**2+(uncrsxscPA*rsrotxscavg/ & rxsPAavg**2)**2)/(1-(rsrotxscavg/rxsPAavg)**2) + & ((unclsxsc/lxsPAavg)**2+(unclsxscPA*lsrotxscavg/ & lxsPAavg**2)**2)/(1-(lsrotxscavg/lxsPAavg)**2) )*1d3 end if c write(3,7006) rsrotgsc/nrgs, c & rsqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs) c & *rSqrt(nrgssq)/nrgs c write(3,7005) lsrotgsc/nlgs c & rsqrt((lsrotgscsq-lsrotgsc**2/nlgs)/nlgs) c & *rSqrt(nlgssq)/nrgs c write(3,7012) rsrotgsc/nrgs - lsrotgsc/nlgs , c & rsqrt( (lsrotgscsq-lsrotgsc**2/nlgs)/nlgs c & *nlgssq/nlgs**2 c & + (rsrotgscsq-rsrotgsc**2/nrgs)/nrgs c & *nrgssq/nrgs**2 ) cccccccc estimated PNC angle scaling TScat and Xscat into average of ALL if (nlxs.gt.0) then rotnsw = dasin(rsrotnscavg/rsnscPAavg)*1d3 urotnsw = rsqrt((uncrsnsc/rsnscPAavg)**2+ & (uncrsnscPA*rsrotnscavg/rsnscPAavg**2)**2)/ & rsqrt(1-(rsrotnscavg/rsnscPAavg)**2)*1d3 rotnse =dasin(lsrotnscavg/lsnscPAavg)*1d3 urotnse = rsqrt((unclsnsc/lsnscPAavg)**2+ !rcm & (unclsnscPA*lsrotnscavg/lsnscPAavg**2)**2)/ & rsqrt(1-(lsrotnscavg/lsnscPAavg)**2)*1d3 rotsw = dasin(rsrottscavg/rtsPAavg)*1d3 urotsw = rsqrt((uncrstsc/rtsPAavg)**2+ & (uncrstscPA*rsrottscavg/rtsPAavg**2)**2)/ & rsqrt(1-(rsrottscavg/rtsPAavg)**2)*1d3 rotse = dasin(lsrottscavg/ltsPAavg)*1d3 urotse = rsqrt((unclstsc/ltsPAavg)**2+ & (unclstscPA*lsrottscavg/ltsPAavg**2)**2)/ & rsqrt(1-(lsrottscavg/ltsPAavg)**2)*1d3 rotxw = dasin(rsrotxscavg/rxsPAavg)*1d3 urotxw = rsqrt((uncrsxsc/rxsPAavg)**2+ & (uncrsxscPA*rsrotxscavg/rxsPAavg**2)**2)/ & rsqrt(1-(rsrotxscavg/rxsPAavg)**2)*1d3 rotxe = dasin(lsrotxscavg/lxsPAavg)*1d3 urotxe = rsqrt((unclsxsc/lxsPAavg)**2+ & (unclsxscPA*lsrotxscavg/lxsPAavg**2)**2)/ & rsqrt(1-(lsrotxscavg/lxsPAavg)**2)*1d3 write(3,*) write(3,*) write(3,889) 0.5*( (fse-fsw)*0.5*(rotnsw+rotnse) + & fsw*rotsw-fse*rotse ), & 0.5*rsqrt( ((0.5*(rotnsw+rotnse)-rotse)*ufse)**2 + & (-(0.5*(rotnsw+rotnse)+rotsw)*ufsw)**2 + & (0.5*(fse-fsw)*urotnse)**2 + (0.5*(fse-fsw)*urotnsw)**2 + & (fsw*urotsw)**2 + (fse*urotse)**2 ) write(3,*) write(3,888) 0.5*( (fse-fsw+fxe-fxw)*0.5*(rotnsw+rotnse) + & fsw*rotsw-fse*rotse + fxw*rotxw-fxe*rotxe ), & 0.5*rsqrt( ((0.5*(rotnsw+rotnse)-rotse)*ufse)**2 + & (-(0.5*(rotnsw+rotnse)+rotsw)*ufsw)**2 + & (0.5*(fse-fsw+fxe-fxw)*urotnse)**2 + & (0.5*(fse-fsw+fxe-fxw)*urotnsw)**2 + & (fsw*urotsw)**2 +(fse*urotse)**2 + & (fxw*urotxw)**2 +(fxe*urotxe)**2 + & (rotxw*ufxw)**2 +(rotxe*ufxe)**2 ) write(3,*) write(3,890) (fse-fsw+fxe-fxw)*0.5*(rotnsw+rotnse), & rsqrt( (ufsw**2+ufse**2+ufxw**2+ufxe**2)* & (0.5*(rotnsw+rotnse))**2 & + ((fse-fsw+fxe-fxw)*0.5)**2*(urotnsw**2+urotnse**2) ) write(3,893) (fse-fsw)*0.5*(rotnsw+rotnse), & rsqrt( (ufsw**2+ufse**2)*(0.5*(rotnsw+rotnse))**2 & + ((fse-fsw)*0.5)**2*(urotnsw**2+urotnse**2)**2 ) write(3,891) fsw*rotsw-fse*rotse , & rsqrt( (ufsw*rotsw)**2 + (fsw*urotsw)**2 & + (ufse*rotse)**2 + (fse*urotse)**2 ) write(3,892) fxw*rotxw-fxe*rotxe , & rsqrt( (ufxw*rotxw)**2 + (fxw*urotxw)**2 & + (ufxe*rotxe)**2 + (fxe*urotxe)**2 ) write(3,*) end if write(3,*) ' ***************************************************', &'******************' write(3,*) if (nder.gt.0) then write(3,702) rsde/nder, & rsqrt((rsdesq-rsde**2/nder)/nder)*rSqrt(ndersq)/nder write(3,701) lsde/ndel, & rsqrt((lsdesq-lsde**2/ndel)/ndel)*rSqrt(ndelsq)/ndel end if if ((nrts.gt.1).and.(nlts.gt.1)) then write(3,802) rwde/nrts, & rsqrt((rwdesq-rwde**2/nrts)/nrts)*rsqrt(nrtssq)/nrts write(3,801) lwde/nlts, & rsqrt((lwdesq-lwde**2/nlts)/nlts)*rsqrt(nltssq)/nlts write(3,822) rpl/nrts, & rsqrt((rplsq-rpl**2/nrts)/nrts)*rsqrt(nrtssq)/nrts write(3,821) lpl/nlts, & rsqrt((lplsq-lpl**2/nlts)/nlts)*rsqrt(nltssq)/nlts end if write(3,*) write(3,*) write(3,*) ' NonScat avg pathlen (cm) West ',rplns/nrnsc,' +- ', & rsqrt((rplnssq-rplns**2/nrnsc)/nrnsc)*rSqrt(nrnscsq)/nrnsc write(3,*) ' NonScat avg pathlen (cm) East ',lplns/nlnsc,' +- ', & rsqrt((lplnssq-lplns**2/nlnsc)/nlnsc)*rSqrt(nlnscsq)/nlnsc write(3,*) write(3,715) nrts, nlts write(3,716) nrts/nlts write(3,718) nrts-nlts,rSqrt(nltssq+nrtssq) write(3,782) nrtsraw,nltsraw write(3,*) write(3,7015) nrgs, nlgs write(3,7016) nrgs/nlgs write(3,7018) nrgs-nlgs,rSqrt(nlgssq+nrgssq) write(3,7082) nrgsraw,nlgsraw write(3,*) write(3,7020) nrxs, nlxs write(3,7021) nrxs/nlxs write(3,7022) nrxs-nlxs,rsqrt(nlxssq+nrxssq) write(3,7023) nrxsraw,nlxsraw write(3,*) if ((nrsc.gt.0).and.(nlsc.gt.0)) then write(3,815) nrsc, nlsc write(3,816) nrsc/nlsc write(3,818) nrsc-nlsc,sqrt(float(nlsc+nrsc)) else write(3,*) ' NO scattering so limited output written ' endif c--------- write(6,*);write(3,*) write(6,*) ' Number that bounced ',numbounce write(3,*) ' Number that bounced ',numbounce write(6,*) ' Number that bounced in 2 guides ',numbounce2 write(3,*) ' Number that bounced in 2 guides ',numbounce2 write(6,*) ' Number that bounced in 3 guides ',numbounce3 write(3,*) ' Number that bounced in 3 guides ',numbounce3 write(6,*) ' Number crossed over after target area WGT: ', & dncross write(3,*) ' Number crossed over after target area WGT: ', & dncross write(3,780) Nxo write(6,780) Nxo write(6,*) write(6,*) ' num guide1 bounced x,y 0x: ',ngbx0,ngby0 write(6,*) ' num guide1 bounced x,y 1x: ',ngbx1,ngby1 write(6,*) ' num guide1 bounced x,y 2x: ',ngbx2,ngby2 write(6,*) ' num guide1 bounced x,y 3x: ',ngbx3,ngby3 write(6,*) ' num guide1 scatters ',gsc write(3,*) write(3,*) ' num guide1 bounced x,y 0x: ',ngbx0,ngby0 write(3,*) ' num guide1 bounced x,y 1x: ',ngbx1,ngby1 write(3,*) ' num guide1 bounced x,y 2x: ',ngbx2,ngby2 write(3,*) ' num guide1 bounced x,y 3x: ',ngbx3,ngby3 write(3,*) ' num guide1 scatters ',gsc write(6,*) write(6,*) ' num incoil bounced x,y 0x: ',nibx0,niby0 write(6,*) ' num incoil bounced x,y 1x: ',nibx1,niby1 write(6,*) ' num incoil bounced x,y 2x: ',nibx2,niby2 write(6,*) ' num incoil bounced x,y 3x: ',nibx3,niby3 write(6,*) ' num incoil scatters ',isc write(3,*) write(3,*) ' num incoil bounced x,y 0x: ',nibx0,niby0 write(3,*) ' num incoil bounced x,y 1x: ',nibx1,niby1 write(3,*) ' num incoil bounced x,y 2x: ',nibx2,niby2 write(3,*) ' num incoil bounced x,y 3x: ',nibx3,niby3 write(3,*) ' num incoil scatters ',isc write(6,*) write(6,*) ' num outcoil bounced x,y 0x: ',nobx0,noby0 write(6,*) ' num outcoil bounced x,y 1x: ',nobx1,noby1 write(6,*) ' num outcoil bounced x,y 2x: ',nobx2,noby2 write(6,*) ' num outcoil bounced x,y 3x: ',nobx3,noby3 write(6,*) ' num outcoil scatters ',osc write(3,*) write(3,*) ' num outcoil bounced x,y 0x: ',nobx0,noby0 write(3,*) ' num outcoil bounced x,y 1x: ',nobx1,noby1 write(3,*) ' num outcoil bounced x,y 2x: ',nobx2,noby2 write(3,*) ' num outcoil bounced x,y 3x: ',nobx3,noby3 write(3,*) ' num outcoil scatters ',osc write(6,*) ; write(3,*) write(6,*) ' Fraction of det hits that bounced ', & numbounce/cot write(3,*) ' Fraction of det hits that bounced ', & numbounce/cot write(6,*) ' Fraction of det hits that bounced in 2 guides', & numbounce2/cot write(3,*) ' Fraction of det hits that bounced in 2 guides ', & numbounce2/cot write(6,*) ' Fraction of det hits that bounced in 3 guides ', & numbounce3/cot write(3,*) ' Fraction of det hits that bounced in 3 guides ', & numbounce3/cot write(6,*) ' Frac. of det hits that crossed after targets WGT ', & dble(dncross)/cot write(3,*) ' Frac. of det hits that crossed after targets WGT', & dble(dncross)/cot write(6,*) ' Frac. of det hits that crossed over anywhere ', & Nxo/cot write(3,*) ' Frac. of det hits that crossed over anywhere ', & Nxo/cot close(unit=3) c bec return end c Functions to avoid domain errors -rcm function rsqrt(B) double precision rsqrt,B nanflg = isnan(B) if(B.lt.0.d0) then write(17,*) ' argsqrt<0 ',B write(17,*) ' argsqrt<0 ',B rsqrt = 0 else if (nanflg.eq..true.) then B = 0 x1 = -1d6; y1 = -1d6 else rsqrt = sqrt(B) end if end function ratan(A) double precision ratan,A,Pi Pi = 3.14159 nanflg = isnan(A) if(A.gt.1.d10) then ratan = Pi/2 else if(A.lt.-1.d10) then ratan = 3/2*Pi else if(nanflg.eq..true.)then ratan = 0 x1 = -1d6; y1 = -1d6 else ratan = atan(A) end if end c function rtan(A) c double precision rtan,A c nanflg = isnan(A) c if(nanflg.eq..true.) then c rtan = 0 c x1 = -1d6; y1 = -1d6 c else c rtan = tan(A) c end if c end c function rsin(A) c double precision rsin,A c nanflg = isnan(A) c if(nanflg.eq..true.) then c rsin = 0 c x1 = -1d6; y1 = -1d6 c else c rsin = sin(A) c end if c end c c function rasin(A) c double precision rasin,A c nanflg = isnan(A) c if(nanflg.eq..true.) then c rasin = 0 c x1 = -1d6; y1 = -1d6 c else c rasin = tan(A) c end if c end c c function rcos(A) c double precision rcos,A c nanflg = isnan(A) c if(nanflg.eq..true.) then c rcos = 0 c x1 = -1d6; y1 = -1d6 c else c rcos = cos(A) c end if c end c c function racos(A) c double precision racos,A c nanflg = isnan(A) c if(nanflg.eq..true.) then c racos = 0 c x1 = -1d6; y1 = -1d6 c else c racos = tan(A) c end if c end c c c