program ffit ! ! Perform fits to magnetic field tables for ! South Arc Upgrade combined-function magnets ! ! Adapted from the 2002 wiggler field fit program ! /home/critten/vf/fieldfit/fitmaind.f ! ! The analytic parameterization is updated to provide ! for arbitrary field shapes. ! ! Internal units are Tesla and meters ! ! Converted from f77 to f90 ! 18 March 2016 ! JAC !-------------------------------------------------------------- implicit none ! real hm common/pawc/hm(800000) ! !------------------------------------------------ integer LUNDAT integer LUNTABLE integer LUNHIS integer LUNTXT common/luncom/LUNDAT,LUNTABLE,LUNHIS,LUNTXT !------------------------------------------------ ! integer icycle,istat !--------------------------------------------------------------- character(*) fname ! parameter ( fname = 'sarc_bq_10_integ_3' ) ! parameter ( fname = 'qf_01_integ' ) parameter ( fname = 'qf_02_integ' ) ! ! Fiducial volume limits (cm) ! --------------------------- real xfid1,yfid1,zfid1,xfid2,yfid2,zfid2 data xfid1,yfid1,zfid1/0.,0.,0./ data xfid2,yfid2,zfid2/3.,1.,15./ ! !--------------------------------------- ! ! Define unit numbers LUNDAT = 30 LUNTABLE = 25 LUNHIS = 60 LUNTXT = 61 ! call hlimit(800000) ! ! Read data from Vector Fields calculation ! Fiducial volume defined here ! fname=fname(1:index(fname,' ')-1) call readfv(xfid1,yfid1,zfid1,xfid2,yfid2,zfid2,fname) ! ! Calculate curl and divergence values call curldiv ! ! Write file of field values within fiducial volume call writefv(fname) ! ! Perform fit call fieldfit ! 99999 continue end subroutine fieldfit ! ! Perform mip spectra fits ! implicit none ! external fcnfield !-------------------------------------------------------- ! Fit parameters ! integer NTERMAX parameter (NTERMAX = 600) integer NP parameter (NP = 7*NTERMAX) character(10) pname(NP) integer nterm double precision par(NP),step(NP),pmin(NP),pmax(NP),sigpar(NP) double precision dbx(np),dby(np),dbz(np) common/fitcom/par,step,pmin,pmax,sigpar,pname,nterm common/dbcom/dbx,dby,dbz !-------------------------------------------------------- ! !------------------------------------------------ integer LUNDAT integer LUNTABLE integer LUNHIS integer LUNTXT common/luncom/LUNDAT,LUNTABLE,LUNHIS,LUNTXT !------------------------------------------------ ! ! double precision fmin,fedm,errdef integer npari,nparx,istat character(80) chdum double precision val,dval,bnd1,bnd2 integer idum ! integer errflag double precision arglis(10) ! real by0x,by0y,by0z external by0x external by0y external by0z ! integer iflag integer i,j,ii,kk ! ! Initialize MINUIT and set unit number print *,' Initializing MINUIT' call mninit(5,6,7) ! call mnseti('Magnetic Field Fits') ! ! Define Parameter names do i=1,ntermax write(pname(7*(i-1)+1),1100)i 1100 format('C',i3.3) write(pname(7*(i-1)+2),1101)i 1101 format('kx',i3.3) write(pname(7*(i-1)+3),1102)i 1102 format('kz',i3.3) write(pname(7*(i-1)+4),1104)i 1104 format('x0',i3.3) write(pname(7*(i-1)+5),1105)i 1105 format('y0',i3.3) write(pname(7*(i-1)+6),1106)i 1106 format('phiz',i3.3) write(pname(7*(i-1)+7),1107)i 1107 format('plane',i3.3) enddo ! ! Initialize parameter array ! Internal units are Tesla and meters call readpar(iflag) ! if(iflag == 1)then write(6,1200) 1200 format(' Routine FIELDFIT: Error from READPAR'/ & ' Using local parameter definitions') ! Use local parameter values appropriate for the QF quadrupole. ! ! Its peak field is about 6 kG. BY is negative at positive X. ! Its transverse good field region is about +-3 cm in both X and Y, ! so kx,ky should be much smaller than pi/0.06 = 50 m-1. ! ! Its longitudinal dependence has a half-period of about 10 cm. ! So kz should be about pi/0.1 = 30 m-1. ! ! It will use primarily the QU family of terms (PLANE=2). nterm=12 do i=1,nterm if(i <= 3)then par(1+(i-1)*7)=-0.6 par(2+(i-1)*7)=1. par(3+(i-1)*7)=30. par(4+(i-1)*7)=0. par(5+(i-1)*7)=0. par(6+(i-1)*7)=0. par(7+(i-1)*7)=2. elseif(i <= 6)then par(1+(i-1)*7)=-0.6 par(2+(i-1)*7)=0.1 par(3+(i-1)*7)=60. par(4+(i-1)*7)=0. par(5+(i-1)*7)=0. par(6+(i-1)*7)=0. par(7+(i-1)*7)=2. elseif(i <= 9)then par(1+(i-1)*7)=-0.6 par(2+(i-1)*7)=1. par(3+(i-1)*7)=60. par(4+(i-1)*7)=0. par(5+(i-1)*7)=0. par(6+(i-1)*7)=0. par(7+(i-1)*7)=2. else par(1+(i-1)*7)=-0.6 par(2+(i-1)*7)=0.1 par(3+(i-1)*7)=30. par(4+(i-1)*7)=0. par(5+(i-1)*7)=0. par(6+(i-1)*7)=0. par(7+(i-1)*7)=2. endif enddo !************ ! Now define step,pmin,pmax for C, kx, kz, x0, y0. do i=1,nterm do j=1,3 step(j+(i-1)*7)=abs(par(j+(i-1)*7)/4.) pmin(j+(i-1)*7)=-5*abs(par(j+(i-1)*7)) pmax(j+(i-1)*7)=5*abs(par(j+(i-1)*7)) enddo ! Fix x0,y0 to zero step(4+(i-1)*7)=0. pmin(4+(i-1)*7)=0. pmax(4+(i-1)*7)=0. step(5+(i-1)*7)=0. pmin(5+(i-1)*7)=0. pmax(5+(i-1)*7)=0. ! Fix PHIZ for maximum transverse field at Z=0 step(6+(i-1)*7)=0. pmin(6+(i-1)*7)=0. pmax(6+(i-1)*7)=0. ! Fix family. Do not vary step(7+(i-1)*7)=0. pmin(7+(i-1)*7)=0. pmax(7+(i-1)*7)=3. enddo endif ! write(6,1250)(par(kk),kk=1,7*nterm) ! 1250 format(' --- Initial Fit Parameters ---'/30(5x,4f14.6/)) 1250 format(' --- Initial Fit Parameters ---'/100(7(2x,e14.7)/)) ! ! Define subset of field points to be used in ! the field calculation call defuse ! ! Write out comparison ntuple call writent ! ! Define Parameters in MINUIT do i=1,nterm*7 ! write (6,1300)i,pname(i),par(i),step(i),pmin(i),pmax(i) 1300 format('Parameter',i3,a10,7f10.3) call mnparm(i,pname(i),par(i),step(i), & pmin(i),pmax(i),errflag) if(errflag /= 0)then write(6,1000)i 1000 format(' !!! Could not define parameter',i3) endif enddo ! ! Set flag to use analytically calculated derivatives ! If the argument=1, MINUIT used these derivatives even if ! they differ from its own numerically calculated values ! Otherwise it uses the analytically calculated values only ! if they agree with its own numerical estimates ! arglis(1)=1. ! When gradients not calculated arglis(1)=0. call mnexcm(fcnfield,'set gradi ',arglis,1,errflag,0) ! ! Set print flag ! arglis(1)=10. arglis(1)=5. call mnexcm(fcnfield,'set print ',arglis,1,errflag,0) ! Set strategy (0: minimize calls, 1: default, 2: many calls) !arglis(1)=0. ! 19 March 2016 arglis(1)=2. call mnexcm(fcnfield,'set str ',arglis,1,errflag,0) ! Perform minimization arglis(1)=20000 ! 20 March 2016 arglis(1)=50000 arglis(2)=0.001 print *,' calling migrad. arglis=',arglis call mnexcm(fcnfield,'migrad ',arglis,2,errflag,0) print *,' done with migrad' ! Get Fit Results call mnstat (fmin,fedm,errdef,npari,nparx,istat) if ( istat < 3)then write(6,8110)istat 8110 format(' MIGRAD did not converge --- istat=',i2) else write(6,8120)istat 8120 format(' MIGRAD converged *** istat=',i2) endif ! ! Print correlations ! arglis(1)=0. ! call mnexcm(fcnfield,'show cor ',arglis,1,errflag,0) ! ! if ( istat == 3 )then do i=1,nterm*7 ! Get parameters and errors call mnpout(i,chdum,val,dval,bnd1,bnd2,idum) par(i)=val sigpar(i)=dval enddo ! else ! Zero parameters for output if no convergence ! do i=1,nterm*7 ! par(i)=0. ! sigpar(i)=0. ! enddo ! endif ! ! Write final parameters to file call writepar(iflag) ! ! Do further error calculation ! call mnexcm(fcnfield,'minos ',0,0,errflag,0) ! ! Write out ntuple file following fit call writent ! call mnexcm(fcnfield,'stop ',arglis,1,errflag,0) ! 999 continue ! end subroutine fcnfield(npar,gin,f,pp,iflag) ! ! MINUIT utility routine for the field fits ! ! Added code for writing out intermediate results ! 28 June 2002 JAC ! implicit none !-------------------------------------------------------- ! Fit parameters ! integer NTERMAX parameter (NTERMAX = 600) integer NP parameter (NP = 7*NTERMAX) character(10) pname(NP) integer nterm double precision par(NP),step(NP),pmin(NP),pmax(NP),sigpar(NP) double precision dbx(np),dby(np),dbz(np) common/fitcom/par,step,pmin,pmax,sigpar,pname,nterm common/dbcom/dbx,dby,dbz !----------------------------------------------------------- ! Arguments ! integer npar double precision f,gin(NP),pp(NP) ! dimension gin(*),pp(*) integer iflag ! !----------------------------------------------------------- ! Common block for fit numbering integer numfit common/fitnum/numfit ! !-------------------------------------------------------- double precision bxfun,byfun,bzfun external bxfun external byfun external bzfun !-------------------------------------------------------- !============================================================ ! Common block of field values from Vector Fields character(80) vfname integer nsize parameter (nsize=1500000) integer nbcalc logical accept(nsize) double precision vfx(nsize) double precision vfbx(nsize) double precision vfebx(nsize) double precision vfy(nsize) double precision vfby(nsize) double precision vfeby(nsize) double precision vfz(nsize) double precision vfbz(nsize) double precision vfebz(nsize) common/vfcom/vfx,vfy,vfz,vfbx,vfby,vfbz, & vfebx,vfeby,vfebz, & nbcalc,accept,vfname !============================================================ ! double precision chi2,scalef double precision x,y,z,bx,by,bz integer nchi2,npp double precision del,bmod2,sumb,relres,wrelres ! real rmsres ! integer i,j,kk,ienter data ienter/0/ ! character(50) fname integer ncall,ncprint data ncall,ncprint/0,200/ ! ienter=ienter+1 ! ! copy parameters to function common block do i=1,np par(i)=pp(i) enddo ! !=================================================== chi2=0. nchi2=0 ! scalef=1.e-7 scalef=1. do i=1,nbcalc ! Get rid of some outlying points if(accept(i))then nchi2=nchi2+1 x=vfx(i) y=vfy(i) z=vfz(i) bx=bxfun(x,y,z) by=byfun(x,y,z) bz=bzfun(x,y,z) chi2=chi2+scalef*((vfbx(i)-bx)/vfebx(i))**2 ! if(i <= 2)print *,' i,x,chi2x=',i,x,chi2 chi2=chi2+scalef*((vfby(i)-by)/vfeby(i))**2 ! if(abs(i-1) <= 0)write(6,6500)i,by,vfby(i),chi2 ! if(abs(i-51) <= 0)write(6,6500)i,by,vfby(i),chi2 6500 format(' i,fby,by,chi2x+chi2y=',i6,2f11.4,f13.6) chi2=chi2+scalef*((vfbz(i)-bz)/vfebz(i))**2 endif enddo ! f=chi2 ! ! Calculate derivatives for IFLAG=2 !----------------------------------- if(iflag == 2)then ! Skip derivative calculation until the functions DBXFUN, ... ! have been updated. NB: Set GRADI should be set to 0 in ! the subroutine FIELDFIT. ! 7 March 2016 jac goto 700 do i=1,7*nterm gin(i)=0. enddo do i=1,nbcalc if(accept(i))then x=vfx(i) y=vfy(i) z=vfz(i) bx=bxfun(x,y,z) by=byfun(x,y,z) bz=bzfun(x,y,z) call dbxfun(x,y,z) call dbyfun(x,y,z) call dbzfun(x,y,z) ! do j=1,7*nterm gin(j)=gin(j) & -dbx(j)*2*scalef*(vfbx(i)-bx)/(vfebx(i)**2) & -dby(j)*2*scalef*(vfby(i)-by)/(vfeby(i)**2) & -dbz(j)*2*scalef*(vfbz(i)-bz)/(vfebz(i)**2) ! if(j == 1.and.i <= 5)write(6,6700)i,x,y,z ! if(j == 1.and.i <= 5)write(6,6710)bx,by,bz ! if(j == 1.and.i <= 5)write(6,6720)vfbx(i),vfby(i),vfbz(i) 6700 format(' C1 derivative for Field point',i6, & ': x,y,z=',3f9.5) 6710 format(' Field function values bx,by,bz=',3f11.4) 6720 format(' Field table values bx,by,bz= ',3f11.4/) enddo endif enddo ! write(6,7000)(kk,gin(kk),kk=1,7*nterm) 7000 format(' Array of derivatives GIN',/100(4(i3,2x,e16.7,2x)/)) 700 continue endif ! ! print *,' ncall,chi2=',ncall+1,f !=================================================== ! ! Write out parameter file with intermediate results ! ! ncall=ncall+1 if(ncall == 1.or.mod(ncall,ncprint) == 0)then write(fname,8000)numfit+1,ncall 8000 format('in/fit',i4.4,'_',i6.6,'.in') open(27,file=fname,err=800,status='new') call parout(27) close(27) ! Write parameters to file in BMAD format write(fname,8010)numfit+1,ncall 8010 format(i4.4,'_',i6.6) fname='bmad/'//vfname(1:index(vfname,' ')-1)//'_'//fname call wbmad(fname) goto 850 800 continue write(6,8100)fname 8100 format(' Routine FCNFIELD could not open file ',a17, & '. Skipped writing intermediate fit result.') 850 continue ! ! Calculate rms residual in gauss rmsres=0. relres=0. wrelres=0. sumb=0. do i=1,nbcalc x=vfx(i) y=vfy(i) z=vfz(i) bx=10000.*bxfun(x,y,z) by=10000.*byfun(x,y,z) bz=10000.*bzfun(x,y,z) del = (10000.*vfbx(i)-bx)**2 & +(10000.*vfby(i)-by)**2 & +(10000.*vfbz(i)-bz)**2 rmsres=rmsres+del bmod2=bx**2+by**2+bz**2 if(bmod2 > 0)then relres=relres+del/bmod2 wrelres=wrelres+del/sqrt(bmod2) sumb=sumb+sqrt(bmod2) endif ! if(mod(i-1,20000) == 0.or.i == 51)then write(6,8500)i,10000.*vfbx(i),10000.*vfbx(i)-bx, & 10000.*vfby(i),10000.*vfby(i)-by, & 10000.*vfbz(i),10000.*vfbz(i)-bz 8500 format(' **FCNFIELD Point',i6,' vfbx,dbx=',2f11.2, & ' vfby,dby=',2f11.2, & ' vfbz,dbz=',2f11.2) endif enddo ! write(6,8050)ncall,sqrt(rmsres),nbcalc,sqrt(rmsres/nbcalc),chi2/scalef 8050 format(' >>FCNFIELD call ',i6.6, & ' Root of sum of squared residuals=',e11.4, & ' gauss for ',i6,' points'/, & ' >>',18x,' RMS residual per point= ', & f15.4,' gauss'/' >>',18x,' Chi-squared=',f13.2,' G^2') ! write(6,8055)100*sqrt(relres)/nbcalc,100*sqrt(wrelres)/sumb 8055 format(' >>',18x,' Relative rms residual per point= ',e10.3,' %' & /,' >>',18x,' Field-weighted relative rms residual per point= ', & e10.3,' %') ! endif ! ! 99999 continue end double precision function bxfun(x,y,z) !********************************************************************** ! ! Fit function for magnetic field ! ! J.A.Crittenden ! 13 March 2002 ! !********************************************************************** implicit none ! ! Input argument double precision x,y,z ! !-------------------------------------------------------- ! Fit parameters ! integer NTERMAX parameter (NTERMAX = 600) integer NP parameter (NP = 7*NTERMAX) character(10) pname(NP) integer nterm double precision par(NP),step(NP),pmin(NP),pmax(NP),sigpar(NP) double precision dbx(np),dby(np),dbz(np) common/fitcom/par,step,pmin,pmax,sigpar,pname,nterm common/dbcom/dbx,dby,dbz !-------------------------------------------------------- integer i,j,ipl double precision term double precision coef,kx,ky,kz,x0,y0,phiz,plane !-------------------------------------------------------- bxfun=0. do i=1,nterm j=7*(i-1) coef=par(j+1) kx=par(j+2) kz=par(j+3) x0=par(j+4) y0=par(j+5) phiz=par(j+6) plane=par(j+7) ipl=nint(plane) if(ipl == 0)then if(kx > 0)then ky=sqrt(kx**2+kz**2) term=(kx/ky)*cos(kx*(x+x0))*cosh(ky*(y+y0))*cos(kz*z+phiz) elseif(abs(kx) < abs(kz))then ky=sqrt(-kx**2+kz**2) term=(kx/kz)*cosh(kx*(x+x0))*cosh(ky*(y+y0))*cos(kz*z+phiz) elseif(abs(kx) >= abs(kz))then ky=sqrt(kx**2-kz**2) term=cosh(kx*(x+x0))*cos(ky*(y+y0))*cos(kz*z+phiz) endif elseif(ipl == 1)then if(kx > 0)then ky=sqrt(kx**2+kz**2) term=-(kx/ky)*sin(kx*(x+x0))*sinh(ky*(y+y0))*cos(kz*z+phiz) elseif(abs(kx) < abs(kz))then ky=sqrt(-kx**2+kz**2) term=(kx/kz)*sinh(kx*(x+x0))*sinh(ky*(y+y0))*cos(kz*z+phiz) elseif(abs(kx) >= abs(kz))then ky=sqrt(kx**2-kz**2) term=sinh(kx*(x+x0))*sin(ky*(y+y0))*cos(kz*z+phiz) endif elseif(ipl == 2)then if(kx > 0)then ky=sqrt(kx**2+kz**2) term=(kx/ky)*cos(kx*(x+x0))*sinh(ky*(y+y0))*cos(kz*z+phiz) elseif(abs(kx) < abs(kz))then ky=sqrt(-kx**2+kz**2) term=(kx/kz)*cosh(kx*(x+x0))*sinh(ky*(y+y0))*cos(kz*z+phiz) elseif(abs(kx) >= abs(kz))then ky=sqrt(kx**2-kz**2) term=cosh(kx*(x+x0))*sin(ky*(y+y0))*cos(kz*z+phiz) endif elseif(ipl == 3)then if(kx > 0)then ky=sqrt(kx**2+kz**2) term=-(kx/ky)*sin(kx*(x+x0))*cosh(ky*(y+y0))*cos(kz*z+phiz) elseif(abs(kx) < abs(kz))then ky=sqrt(-kx**2+kz**2) term=(kx/kz)*sinh(kx*(x+x0))*cosh(ky*(y+y0))*cos(kz*z+phiz) elseif(abs(kx) >= abs(kz))then ky=sqrt(kx**2-kz**2) term=-sinh(kx*(x+x0))*cos(ky*(y+y0))*cos(kz*z+phiz) endif else print *,' --- Illegal plane value ',ipl,' Stopping.' stop endif bxfun=bxfun+coef*term enddo ! end double precision function byfun(x,y,z) !********************************************************************** ! ! Fit function for magnetic field ! ! J.A.Crittenden ! 13 March 2002 ! !********************************************************************** implicit none ! ! Input argument double precision x,y,z ! !-------------------------------------------------------- ! Fit parameters ! integer NTERMAX parameter (NTERMAX = 600) integer NP parameter (NP = 7*NTERMAX) character(10) pname(NP) integer nterm double precision par(NP),step(NP),pmin(NP),pmax(NP),sigpar(NP) double precision dbx(np),dby(np),dbz(np) common/fitcom/par,step,pmin,pmax,sigpar,pname,nterm common/dbcom/dbx,dby,dbz !-------------------------------------------------------- integer i,j,ipl double precision term double precision coef,kx,ky,kz,x0,y0,phiz,plane !-------------------------------------------------------- byfun=0. do i=1,nterm j=7*(i-1) coef=par(j+1) kx=par(j+2) kz=par(j+3) x0=par(j+4) y0=par(j+5) phiz=par(j+6) plane=par(j+7) ipl=nint(plane) if(ipl == 0)then if(kx > 0)then ky=sqrt(kx**2+kz**2) term=sin(kx*(x+x0))*sinh(ky*(y+y0))*cos(kz*z+phiz) elseif(abs(kx) < abs(kz))then ky=sqrt(-kx**2+kz**2) term=(ky/kz)*sinh(kx*(x+x0))*sinh(ky*(y+y0))*cos(kz*z+phiz) elseif(abs(kx) >= abs(kz))then ky=sqrt(kx**2-kz**2) term=-(ky/kx)*sinh(kx*(x+x0))*sin(ky*(y+y0))*cos(kz*z+phiz) endif elseif(ipl == 1)then if(kx > 0)then ky=sqrt(kx**2+kz**2) term=cos(kx*(x+x0))*cosh(ky*(y+y0))*cos(kz*z+phiz) elseif(abs(kx) < abs(kz))then ky=sqrt(-kx**2+kz**2) term=(ky/kz)*cosh(kx*(x+x0))*cosh(ky*(y+y0))*cos(kz*z+phiz) elseif(abs(kx) >= abs(kz))then ky=sqrt(kx**2-kz**2) term=(ky/kx)*cosh(kx*(x+x0))*cos(ky*(y+y0))*cos(kz*z+phiz) endif elseif(ipl == 2)then if(kx > 0)then ky=sqrt(kx**2+kz**2) term=sin(kx*(x+x0))*cosh(ky*(y+y0))*cos(kz*z+phiz) elseif(abs(kx) < abs(kz))then ky=sqrt(-kx**2+kz**2) term=(ky/kz)*sinh(kx*(x+x0))*cosh(ky*(y+y0))*cos(kz*z+phiz) elseif(abs(kx) >= abs(kz))then ky=sqrt(kx**2-kz**2) term=(ky/kx)*sinh(kx*(x+x0))*cos(ky*(y+y0))*cos(kz*z+phiz) endif elseif(ipl == 3)then if(kx > 0)then ky=sqrt(kx**2+kz**2) term=cos(kx*(x+x0))*sinh(ky*(y+y0))*cos(kz*z+phiz) elseif(abs(kx) < abs(kz))then ky=sqrt(-kx**2+kz**2) term=(ky/kz)*cosh(kx*(x+x0))*sinh(ky*(y+y0))*cos(kz*z+phiz) elseif(abs(kx) >= abs(kz))then ky=sqrt(kx**2-kz**2) term=(ky/kx)*cosh(kx*(x+x0))*sin(ky*(y+y0))*cos(kz*z+phiz) endif else print *,' --- Illegal plane value ',ipl,' Stopping.' stop endif byfun=byfun+coef*term enddo end double precision function bzfun(x,y,z) !********************************************************************** ! ! Fit function for magnetic field ! ! J.A.Crittenden ! 13 March 2002 ! !********************************************************************** implicit none ! ! Input argument double precision x,y,z ! !-------------------------------------------------------- ! Fit parameters ! integer NTERMAX parameter (NTERMAX = 600) integer NP parameter (NP = 7*NTERMAX) character(10) pname(NP) integer nterm double precision par(NP),step(NP),pmin(NP),pmax(NP),sigpar(NP) double precision dbx(np),dby(np),dbz(np) common/fitcom/par,step,pmin,pmax,sigpar,pname,nterm common/dbcom/dbx,dby,dbz !-------------------------------------------------------- integer i,j,ipl double precision term double precision coef,kx,ky,kz,x0,y0,phiz,plane !-------------------------------------------------------- bzfun=0. do i=1,nterm j=7*(i-1) coef=par(j+1) kx=par(j+2) kz=par(j+3) x0=par(j+4) y0=par(j+5) phiz=par(j+6) plane=par(j+7) ipl=nint(plane) if(ipl == 0)then if(kx > 0)then ky=sqrt(kx**2+kz**2) term=-(kz/ky)*sin(kx*(x+x0))*cosh(ky*(y+y0))*sin(kz*z+phiz) elseif(abs(kx) < abs(kz))then ky=sqrt(-kx**2+kz**2) term=-sinh(kx*(x+x0))*cosh(ky*(y+y0))*sin(kz*z+phiz) elseif(abs(kx) >= abs(kz))then ky=sqrt(kx**2-kz**2) term=-(kz/kx)*sinh(kx*(x+x0))*cos(ky*(y+y0))*sin(kz*z+phiz) endif elseif(ipl == 1)then if(kx > 0)then ky=sqrt(kx**2+kz**2) term=-(kz/ky)*cos(kx*(x+x0))*sinh(ky*(y+y0))*sin(kz*z+phiz) elseif(abs(kx) < abs(kz))then ky=sqrt(-kx**2+kz**2) term=-cosh(kx*(x+x0))*sinh(ky*(y+y0))*sin(kz*z+phiz) elseif(abs(kx) >= abs(kz))then ky=sqrt(kx**2-kz**2) term=-(kz/kx)*cosh(kx*(x+x0))*sin(ky*(y+y0))*sin(kz*z+phiz) endif elseif(ipl == 2)then if(kx > 0)then ky=sqrt(kx**2+kz**2) term=-(kz/ky)*sin(kx*(x+x0))*sinh(ky*(y+y0))*sin(kz*z+phiz) elseif(abs(kx) < abs(kz))then ky=sqrt(-kx**2+kz**2) term=-sinh(kx*(x+x0))*sinh(ky*(y+y0))*sin(kz*z+phiz) elseif(abs(kx) >= abs(kz))then ky=sqrt(kx**2-kz**2) term=-(kz/kx)*sinh(kx*(x+x0))*sin(ky*(y+y0))*sin(kz*z+phiz) endif elseif(ipl == 3)then if(kx > 0)then ky=sqrt(kx**2+kz**2) term=-(kz/ky)*cos(kx*(x+x0))*cosh(ky*(y+y0))*sin(kz*z+phiz) elseif(abs(kx) < abs(kz))then ky=sqrt(-kx**2+kz**2) term=-cosh(kx*(x+x0))*cosh(ky*(y+y0))*sin(kz*z+phiz) elseif(abs(kx) >= abs(kz))then ky=sqrt(kx**2-kz**2) term=-(kz/kx)*cosh(kx*(x+x0))*cos(ky*(y+y0))*sin(kz*z+phiz) endif else print *,' --- Illegal plane value ',ipl,' Stopping.' stop endif bzfun=bzfun+coef*term enddo ! end subroutine dbxfun(x,y,z) !********************************************************************** ! ! Calculates derivatives of fit function for each parameter and ! stores them in common block DBCOM for use by FCNFIELD ! ! J.A.Crittenden ! 18 April 2002 ! !********************************************************************** implicit none ! ! Input argument double precision x,y,z ! !-------------------------------------------------------- ! Fit parameters ! integer NTERMAX parameter (NTERMAX = 600) integer NP parameter (NP = 7*NTERMAX) character(10) pname(NP) integer nterm double precision par(NP),step(NP),pmin(NP),pmax(NP),sigpar(NP) double precision dbx(np),dby(np),dbz(np) common/fitcom/par,step,pmin,pmax,sigpar,pname,nterm common/dbcom/dbx,dby,dbz !-------------------------------------------------------- integer i,j double precision term,dterm double precision coef,kx,ky,kz,x0,y0,phiz,plane !-------------------------------------------------------- do i=1,nterm j=7*(i-1) coef=par(j+1) kx=par(j+2) kz=par(j+3) x0=par(j+4) y0=par(j+5) phiz=par(j+6) plane=par(j+7) ! if(kx > 0)then ky=sqrt(kx**2+kz**2) term=sin(kx*(x+x0))*sinh(ky*(y+y0))*cos(kz*z+phiz) ! calculate the four derivatives for kx>0 ! dbx/dcoef dbx(j+1)=(-kx/ky)*term ! dbx/dkx dterm=x*cos(kx*(x+x0))*sinh(ky*(y+y0))*cos(kz*z+phiz) dterm=dterm+sin(kx*(x+x0))*(y+y0)*(kx/ky)*cosh(ky*(y+y0))*cos(kz*z+phiz) dterm=(-kx*dterm-(1-(kx/ky)**2)*term)/ky dbx(j+2)=coef*dterm ! dbx/dkz dterm=sin(kx*(x+x0))*sinh(ky*(y+y0))*(-z*sin(kz*z+phiz)) dterm=dterm+sin(kx*(x+x0))*(y+y0)*(kz/ky)*cosh(ky*(y+y0))*cos(kz*z+phiz) dterm=(-kx*dterm+(kx*kz/(ky**2))*term)/ky dbx(j+3)=coef*dterm ! dbx/dphiz dterm=-sin(kx*(x+x0))*sinh(ky*(y+y0))*sin(kz*z+phiz) dbx(j+4)=coef*(-kx/ky)*dterm elseif(abs(kx) < abs(kz))then ky=sqrt(-kx**2+kz**2) term=-sinh(kx*(x+x0))*sinh(ky*(y+y0))*cos(kz*z+phiz) ! calculate the four derivatives for abs(kx)= abs(kz))then ky=sqrt(kx**2-kz**2) term=-sinh(kx*(x+x0))*sin(ky*(y+y0))*cos(kz*z+phiz) ! calculate the four derivatives for abs(kx)>=abs(kz) ! dbx/dcoef dbx(j+1)=(-kx/ky)*term ! dbx/dkx dterm=-x*cosh(kx*(x+x0))*sin(ky*(y+y0))*cos(kz*z+phiz) dterm=dterm-sinh(kx*(x+x0))*(y+y0)*(kx/ky)*cos(ky*(y+y0))*cos(kz*z+phiz) dterm=(-kx*dterm-(1-(kx/ky)**2)*term)/ky dbx(j+2)=coef*dterm ! dbx/dkz dterm=-sinh(kx*(x+x0))*sin(ky*(y+y0))*(-z*sin(kz*z+phiz)) dterm=dterm-sinh(kx*(x+x0))*(y+y0)*(-kz/ky)*cos(ky*(y+y0))*cos(kz*z+phiz) dterm=(-kx*dterm-(kx*kz/(ky**2))*term)/ky dbx(j+3)=coef*dterm ! dbx/dphiz dterm=sinh(kx*(x+x0))*sin(ky*(y+y0))*sin(kz*z+phiz) dbx(j+4)=coef*(-kx/ky)*dterm endif ! enddo ! end subroutine dbyfun(x,y,z) !********************************************************************** ! ! Calculates derivatives of fit function for each parameter and ! stores them in common block DBCOM for use by FCNFIELD ! ! J.A.Crittenden ! 18 April 2002 ! !********************************************************************** implicit none ! ! Input argument double precision x,y,z ! !-------------------------------------------------------- ! Fit parameters ! integer NTERMAX parameter (NTERMAX = 600) integer NP parameter (NP = 7*NTERMAX) character(10) pname(NP) integer nterm double precision par(NP),step(NP),pmin(NP),pmax(NP),sigpar(NP) double precision dbx(np),dby(np),dbz(np) common/fitcom/par,step,pmin,pmax,sigpar,pname,nterm common/dbcom/dbx,dby,dbz !-------------------------------------------------------- integer i,j double precision term,dterm double precision coef,kx,ky,kz,x0,y0,phiz,plane !-------------------------------------------------------- do i=1,nterm j=7*(i-1) coef=par(j+1) kx=par(j+2) kz=par(j+3) x0=par(j+4) y0=par(j+5) phiz=par(j+6) plane=par(j+7) ! if(kx > 0)then ky=sqrt(kx**2+kz**2) term=cos(kx*(x+x0))*cosh(ky*(y+y0))*cos(kz*z+phiz) ! calculate four derivatives for kx>0 ! dby/dcoef dby(j+1)=term ! dby/dkx dterm=-x*sin(kx*(x+x0))*cosh(ky*(y+y0))*cos(kz*z+phiz) dterm=dterm+cos(kx*(x+x0))*(y+y0)*(kx/ky)*sinh(ky*(y+y0))*cos(kz*z+phiz) dby(j+2)=coef*dterm ! dby/dkz dterm=cos(kx*(x+x0))*cosh(ky*(y+y0))*(-z*sin(kz*z+phiz)) dterm=dterm+cos(kx*(x+x0))*(y+y0)*(kz/ky)*sinh(ky*(y+y0))*cos(kz*z+phiz) dby(j+3)=coef*dterm ! dby/dphiz dterm=-cos(kx*(x+x0))*cosh(ky*(y+y0))*sin(kz*z+phiz) dby(j+4)=coef*dterm elseif(abs(kx) < abs(kz))then ky=sqrt(-kx**2+kz**2) term=cosh(kx*(x+x0))*cosh(ky*(y+y0))*cos(kz*z+phiz) ! calculate four derivatives for abs(kx)= abs(kz))then ky=sqrt(kx**2-kz**2) term=cosh(kx*(x+x0))*cos(ky*(y+y0))*cos(kz*z+phiz) ! calculate four derivatives for abs(kx)>=abs(kz) ! dby/dcoef dby(j+1)=term ! dby/dkx dterm=x*sinh(kx*(x+x0))*cos(ky*(y+y0))*cos(kz*z+phiz) dterm=dterm-cosh(kx*(x+x0))*(y+y0)*(kx/ky)*sin(ky*(y+y0))*cos(kz*z+phiz) dby(j+2)=coef*dterm ! dby/dkz dterm=cosh(kx*(x+x0))*cos(ky*(y+y0))*(-z*sin(kz*z+phiz)) dterm=dterm-cosh(kx*(x+x0))*(y+y0)*(-kz/ky)*sin(ky*(y+y0))*cos(kz*z+phiz) dby(j+3)=coef*dterm ! dby/dphiz dterm=-cosh(kx*(x+x0))*cos(ky*(y+y0))*sin(kz*z+phiz) dby(j+4)=coef*dterm endif ! enddo ! end subroutine dbzfun(x,y,z) !********************************************************************** ! ! Calculates derivatives for each parameter as well and ! stores them in common block DBCOM for use by FCNFIELD ! ! J.A.Crittenden ! 18 April 2002 ! !********************************************************************** implicit none ! ! Input argument double precision x,y,z ! !-------------------------------------------------------- ! Fit parameters ! integer NTERMAX parameter (NTERMAX = 600) integer NP parameter (NP = 7*NTERMAX) character(10) pname(NP) integer nterm double precision par(NP),step(NP),pmin(NP),pmax(NP),sigpar(NP) double precision dbx(np),dby(np),dbz(np) common/fitcom/par,step,pmin,pmax,sigpar,pname,nterm common/dbcom/dbx,dby,dbz !-------------------------------------------------------- integer i,j double precision term,dterm double precision coef,kx,ky,kz,x0,y0,phiz,plane !-------------------------------------------------------- do i=1,nterm j=7*(i-1) j=7*(i-1) coef=par(j+1) kx=par(j+2) kz=par(j+3) x0=par(j+4) y0=par(j+5) phiz=par(j+6) plane=par(j+7) ! if(kx > 0)then ky=sqrt(kx**2+kz**2) term=cos(kx*(x+x0))*sinh(ky*(y+y0))*sin(kz*z+phiz) ! calculate four derivatives for kx>0 ! dbz/dcoef dbz(j+1)=(-kz/ky)*term ! dbz/dkx dterm=-x*sin(kx*(x+x0))*sinh(ky*(y+y0))*sin(kz*z+phiz) dterm=dterm+cos(kx*(x+x0))*(y+y0)*(kx/ky)*cosh(ky*(y+y0))*sin(kz*z+phiz) dterm=(-kz*dterm+(kx*kz/(ky**2))*term)/ky dbz(j+2)=coef*dterm ! dbz/dkz dterm=cos(kx*(x+x0))*sinh(ky*(y+y0))*(z*cos(kz*z+phiz)) dterm=dterm+cos(kx*(x+x0))*(y+y0)*(kz/ky)*cosh(ky*(y+y0))*sin(kz*z+phiz) dterm=(-kz*dterm-(1-(kz/ky)**2)*term)/ky dbz(j+3)=coef*dterm ! dbz/dphiz dterm=cos(kx*(x+x0))*sinh(ky*(y+y0))*cos(kz*z+phiz) dbz(j+4)=coef*(-kz/ky)*dterm elseif(abs(kx) < abs(kz))then ky=sqrt(-kx**2+kz**2) term=cosh(kx*(x+x0))*sinh(ky*(y+y0))*sin(kz*z+phiz) ! calculate four derivatives for abs(kx)= abs(kz))then ky=sqrt(kx**2-kz**2) term=cosh(kx*(x+x0))*sin(ky*(y+y0))*sin(kz*z+phiz) ! calculate dterm for abs(kx)>=abs(kz) ! dbz/dcoef dbz(j+1)=(-kz/ky)*term ! dbz/dkx dterm=x*sinh(kx*(x+x0))*sin(ky*(y+y0))*sin(kz*z+phiz) dterm=dterm+cosh(kx*(x+x0))*(y+y0)*(kx/ky)*cos(ky*(y+y0))*sin(kz*z+phiz) ! sign changed here 6 may 2002 dterm=(-kz*dterm+(kx*kz/(ky**2))*term)/ky dbz(j+2)=coef*dterm ! dbz/dkz dterm=cosh(kx*(x+x0))*sin(ky*(y+y0))*(z*cos(kz*z+phiz)) dterm=dterm+cosh(kx*(x+x0))*(y+y0)*(-kz/ky)*cos(ky*(y+y0))*sin(kz*z+phiz) dterm=(-kz*dterm-(1+(kz/ky)**2)*term)/ky dbz(j+3)=coef*dterm ! dbz/dphiz dterm=cosh(kx*(x+x0))*sin(ky*(y+y0))*cos(kz*z+phiz) dbz(j+4)=coef*(-kz/ky)*dterm endif ! enddo ! end subroutine fillnt ! ! Book and fill ntuple with fit results ! ! J.A.Crittenden ! March 2002 ! implicit none !-------------------------------------------------------- ! Fit parameters ! integer NTERMAX parameter (NTERMAX = 600) integer NP parameter (NP = 7*NTERMAX) character(10) pname(NP) integer nterm double precision par(NP),step(NP),pmin(NP),pmax(NP),sigpar(NP) double precision dbx(np),dby(np),dbz(np) common/fitcom/par,step,pmin,pmax,sigpar,pname,nterm common/dbcom/dbx,dby,dbz !-------------------------------------------------------- integer LUNDAT integer LUNTABLE integer LUNHIS integer LUNTXT common/luncom/LUNDAT,LUNTABLE,LUNHIS,LUNTXT !------------------------------------------------ !============================================================ ! Common block of field values from Vector Fields character(80) vfname integer nsize parameter (nsize=1500000) integer nbcalc logical accept(nsize) double precision vfx(nsize) double precision vfbx(nsize) double precision vfebx(nsize) double precision vfy(nsize) double precision vfby(nsize) double precision vfeby(nsize) double precision vfz(nsize) double precision vfbz(nsize) double precision vfebz(nsize) common/vfcom/vfx,vfy,vfz,vfbx,vfby,vfbz, & vfebx,vfeby,vfebz, & nbcalc,accept,vfname !============================================================ !============================================================ ! Common block of curl and divergence values double precision div(nsize) double precision curlx(nsize) double precision curly(nsize) double precision curlz(nsize) common/curlcom/div,curlx,curly,curlz !============================================================ ! double precision bxfun,byfun,bzfun external bxfun,byfun,bzfun ! logical hexist external hexist ! integer i,id0,nvar,icycle character(12) htitle character(8) tags(16) real fn(16) ! data nvar/16/ data htitle/' Field Fit '/ data tags/' x ',' y ',' z ', & ' vfbx ',' vfby ',' vfbz ', & ' vfebx ',' vfeby ',' vfebz ', & ' fbx ',' fby ',' fbz ', & ' div ',' curlx',' curly',' curlz'/ ! ! Booking ! id0=10 if(hexist(id0))then print *,' Rtn FILLNT deleting pre-existing ntuple ',id0 call hdelet(id0) endif call hbookn(id0,htitle,nvar,'//NT',10000,tags) ! ! Filling ! ! Use units of gauss and cm in the ntuple do i=1,nbcalc ! fn(1)=vfx(i)*100. fn(2)=vfy(i)*100. fn(3)=vfz(i)*100. fn(4)=vfbx(i)*10000. fn(5)=vfby(i)*10000. fn(6)=vfbz(i)*10000. fn(7)=vfebx(i)*10000. fn(8)=vfeby(i)*10000. fn(9)=vfebz(i)*10000. fn(10)=bxfun(vfx(i),vfy(i),vfz(i))*10000. fn(11)=byfun(vfx(i),vfy(i),vfz(i))*10000. fn(12)=bzfun(vfx(i),vfy(i),vfz(i))*10000. fn(13)=div(i) fn(14)=curlx(i) fn(15)=curly(i) fn(16)=curlz(i) ! ! if(accept(i))call hfn(id0,fn) call hfn(id0,fn) ! enddo ! print *,' Rtn FILLNT finished filling ntuple ID=',id0 call hprint(id0) ! 99999 continue end SUBROUTINE ROMIN8(FUN,XMIN,XMAX,XINT,DXINT,NITER) !*** !*** ROUTINE INTEGRATES ONE-VARIABLE FUNCTION FUN BETWEEN !*** XMIN AND XMAX. FUN MUST BE DEFINED IN THE WHOLE INTERVAL, !*** ALSO AT THE LIMITS. ALGORITHM USED IS ROMBERG. ALL !*** ARGUMENTS OF ROMIN8 AND CALCULATIONS IN REAL*8. MAXIMUM !*** NUMBER OF FCN. CALLS IS 2**(NITER+1)+1 (DEFAULT = 1025). !*** IN : FUN = FUNCTION TO BE INTEGRATED. MUST BE DECLARED !*** EXTERNAL IN CALLINHG SEQUENCE !*** XMIN = LOWER BOUNDARY OF INTEGRATION INTERVAL !*** XMAX = UPPER BOUNDARY OF INTEGRATION INTERVAL !*** DXINT = DESIRED RELATIVE ACCURACY OF INTEGRATION !*** DXINT = 0 MEANS HIGHEST POSSIBLE ACCURACY !*** NITER = MAX. NUMBER OF ITERATION LOOPS, DEFAULT = 9 !*** INTERNALLY NITER IS LIMITED TO NITER <= 15 !*** OUT: XINT = INTEGRATION RESULT !*** DXINT = ESTIMATE OF ABSOLUTE INTEGRATION ERROR !*** IMPLICIT REAL*8 (A-H,O-Z) REAL*8 SUM(15) DATA NDEF/9/ !----------------------------------------------------------------------- ! CHECK INTEGRATION INTERVAL, SET START VALUES !----------------------------------------------------------------------- CALL NOARG(NARG) IF(NARG <= 5) THEN NITE=NDEF ELSE NITE=MIN0(NITER,15) END IF XINT=0.D0 XMI=XMIN XMA=XMAX IF(XMI == XMA) RETURN SIGN=1.D0 IF(XMI > XMA) THEN HELP=XMI XMI=XMA XMA=HELP SIGN=-1.D0 END IF NBIN=2 DX=(XMA-XMI)/2.D0 XM=XMI XP=XMA FSUM=-(FUN(XM)+FUN(XP))/2.D0 !----------------------------------------------------------------------- ! MAIN LOOP: VARY BIN WIDTHS, ... !----------------------------------------------------------------------- IFLG=0 DO 10 ISTEP=1,NITE DX=DX/2.D0 NBIN=NBIN*2 !----------------------------------------------------------------------- ! ..., CALCULATE FUNCTION VALUES, ... !----------------------------------------------------------------------- X=XMI-DX NDIF=2**(NITE-ISTEP) DO 20 I=0,NBIN X=X+DX IFLG=-IFLG IF(IFLG < 0) GO TO 20 KBIN=1+I*NDIF ARG=X FSUM=FSUM+FUN(ARG) 20 CONTINUE IFLG=1 !----------------------------------------------------------------------- ! ..., CALCULATE TRAPEZE SUMS, ... !----------------------------------------------------------------------- SUM(ISTEP)=FSUM*DX !----------------------------------------------------------------------- ! ..., EXTRAPOLATE TO BIN WIDTH = 0, ... !----------------------------------------------------------------------- K=1 DO 40 I=ISTEP-1,1,-1 K=K*4 40 SUM(I)=SUM(I+1)+(SUM(I+1)-SUM(I))/DBLE(K-1) !----------------------------------------------------------------------- ! ..., AND ESTIMATE ERROR AND DECIDE ABOUT ITERATION STOP !----------------------------------------------------------------------- IF(SUM(1) /= 0.) THEN ERR=DABS((XINT-SIGN*SUM(1))/SUM(1))/4. ELSE ERR=1.D0 END IF IF(ERR < DXINT) THEN XINT=SIGN*SUM(1) GO TO 100 END IF 10 XINT=SIGN*SUM(1) 100 DXINT=ERR*XINT RETURN END !======================================================================= !*** SUBROUTINE ROMINT(FUN,XMIN,XMAX,XINT,DXINT,NITER) !*** !*** ROUTINE INTEGRATES ONE-VARIABLE FUNCTION FUN BETWEEN !*** XMIN AND XMAX. FUN MUST BE DEFINED IN THE WHOLE INTERVAL, !*** ALSO AT THE LIMITS. ALGORITHM USED IS ROMBERG. MAXIMUM !*** NUMBER OF FCN. CALLS IS 2**(NITER+1)+1 (DEFAULT = 1025). !*** IN : FUN = FUNCTION TO BE INTEGRATED. MUST BE DECLARED !*** EXTERNAL IN CALLINHG SEQUENCE !*** XMIN = LOWER BOUNDARY OF INTEGRATION INTERVAL !*** XMAX = UPPER BOUNDARY OF INTEGRATION INTERVAL !*** DXINT = DESIRED RELATIVE ACCURACY OF INTEGRATION !*** DXINT = 0 MEANS HIGHEST POSSIBLE ACCURACY !*** NITER = MAX. NUMBER OF ITERATION LOOPS, DEFAULT = 9 !*** INTERNALLY NITER IS LIMITED TO NITER <= 15 !*** OUT: XINT = INTEGRATION RESULT !*** DXINT = ESTIMATE OF ABSOLUTE INTEGRATION ERROR !*** REAL*8 XMI,XMA,DX,X,HELP,FSUM,SUM(15) DATA NDEF/9/ !----------------------------------------------------------------------- ! CHECK INTEGRATION INTERVAL, SET START VALUES !----------------------------------------------------------------------- CALL NOARG(NARG) IF(NARG <= 5) THEN NITE=NDEF ELSE NITE=MIN0(NITER,15) END IF XINT=0. XMI=XMIN XMA=XMAX IF(XMI == XMA) RETURN SIGN=1. IF(XMI > XMA) THEN HELP=XMI XMI=XMA XMA=HELP SIGN=-1. END IF NBIN=2 DX=(XMA-XMI)/2. XM=XMI XP=XMA FSUM=-(FUN(XM)+FUN(XP))/2. !----------------------------------------------------------------------- ! MAIN LOOP: VARY BIN WIDTHS, ... !----------------------------------------------------------------------- IFLG=0 DO 10 ISTEP=1,NITE DX=DX/2. NBIN=NBIN*2 !----------------------------------------------------------------------- ! ..., CALCULATE FUNCTION VALUES, ... !----------------------------------------------------------------------- X=XMI-DX NDIF=2**(NITE-ISTEP) DO 20 I=0,NBIN X=X+DX IFLG=-IFLG IF(IFLG < 0) GO TO 20 KBIN=1+I*NDIF ARG=X FSUM=FSUM+FUN(ARG) 20 CONTINUE IFLG=1 !----------------------------------------------------------------------- ! ..., CALCULATE TRAPEZE SUMS, ... !----------------------------------------------------------------------- SUM(ISTEP)=FSUM*DX !----------------------------------------------------------------------- ! ..., EXTRAPOLATE TO BIN WIDTH = 0, ... !----------------------------------------------------------------------- K=1 DO 40 I=ISTEP-1,1,-1 K=K*4 40 SUM(I)=SUM(I+1)+(SUM(I+1)-SUM(I))/DBLE(K-1) !----------------------------------------------------------------------- ! ..., AND ESTIMATE ERROR AND DECIDE ABOUT ITERATION STOP !----------------------------------------------------------------------- IF(SUM(1) /= 0.) THEN ERR=DABS((XINT-SIGN*SUM(1))/SUM(1))/4. ELSE ERR=1. END IF IF(ERR < DXINT) THEN XINT=SIGN*SUM(1) GO TO 100 END IF 10 XINT=SIGN*SUM(1) 100 DXINT=ERR*XINT RETURN END !======================================================================= subroutine hfillf(id,func) ! ! Fill histogram id with values from function func ! implicit none ! integer id real func external func ! character(80) chtit integer nx,ny,nwi,loc real xmi,xma,ymi,yma real dx,x,y integer i ! real con(1000),err(1000) ! logical hexist external hexist ! if (.not. hexist(id)) goto 900 ! call hgive (id,chtit,nx,xmi,xma,ny,ymi,yma,nwi,loc) ! dx = (xma - xmi) / nx ! do i = 1,nx x = xmi + (i-0.5)*dx y = func (x) ! call hf1(id,x,y) con(i)=y err(i)=0. enddo call hpak(id,con) call hpake(id,err) ! goto 99999 900 continue write(6,9000)id 9000 format(' === Routine HFILLF: Histo ID',i5,' does not exist ===') ! 99999 continue end subroutine readfv(xfid1,yfid1,zfid1,xfid2,yfid2,zfid2,fname) !* !* Read field values from file !* created by Vector Fields !* !**************************************************** !* Vector Fields Table File Format !* -------------------------------- !* Units are cm and gauss !* ---------------------- !* 1) The first record contains the number of !* bins in x, y, and z !* 2) The next seven records are strings specifying !* the names of the variables and the units !* 3) The following nx*ny*nz records each !* contain ten free-format numbers, which are !* x, y, z, bx, by, bz, errbx, errby, errbz, mu !* ----------- !* 11 Feb 2002 !* JAC !* ----------- !* Modified to restrict data to fidicial volume !* ---------- !* 13 March 2002 !**************************************************** !* !* Fiducial volume defintion real xfid1,yfid1,zfid1,xfid2,yfid2,zfid2 character(*) fname !* !============================================================ ! Common block of field values from Vector Fields character(80) vfname integer nsize parameter (nsize=1500000) integer nbcalc logical accept(nsize) double precision vfx(nsize) double precision vfbx(nsize) double precision vfebx(nsize) double precision vfy(nsize) double precision vfby(nsize) double precision vfeby(nsize) double precision vfz(nsize) double precision vfbz(nsize) double precision vfebz(nsize) common/vfcom/vfx,vfy,vfz,vfbx,vfby,vfbz, & vfebx,vfeby,vfebz, & nbcalc,accept,vfname !============================================================ character(40) fn character(80) string real vec(10) integer iuse ! Load table file name prefix into common block vfname=fname fn=fname//'.table' print *,' Routine READFV opening file ',fn open(25,file=fn,status='old') read(25,*)nx,ny,nz 1000 format(3i6) ! print *,nx,ny,nz ! do i=1,11 do i=1,7 read(25,2000)string 2000 format(a80) print *,string enddo nbcalc=0 iuse=0 irec=0 10 continue irec=irec+1 ! read(25,*,end=900)(vec(kk),kk=1,10) read(25,*,end=900)(vec(kk),kk=1,6) ! if(mod(irec-1,1000) == 0)print *,'Record',irec,vec if( & abs(vec(1)) >= xfid1.and. & abs(vec(2)) >= yfid1.and. & abs(vec(3)) >= zfid1.and. & abs(vec(1)) <= xfid2.and. & abs(vec(2)) <= yfid2.and. & abs(vec(3)) <= zfid2 & )then nbcalc=nbcalc+1 vfx(nbcalc)=vec(1)/100. vfy(nbcalc)=vec(2)/100. vfz(nbcalc)=vec(3)/100. vfbx(nbcalc)=vec(4)/10000. vfby(nbcalc)=vec(5)/10000. vfbz(nbcalc)=vec(6)/10000. ! vfebx(nbcalc)=vec(7)/10000. ! vfeby(nbcalc)=vec(8)/10000. ! vfebz(nbcalc)=vec(9)/10000. vfebx(nbcalc)=1. vfeby(nbcalc)=1. vfebz(nbcalc)=1. endif goto 10 900 continue close(25) nprod=nx*ny*nz write(6,8000)irec-1,nbcalc,xfid1,yfid1,zfid1,xfid2,yfid2,zfid2 8000 format(' Routine READFV found',i7,' records and stored', & i7,' field measurements within the fiducial volume.'/ & ' The fiducial volume is bounded by x,y,z greater than:',3f7.2,/ & ' and less than:',3f7.2) ! end subroutine writefv(fname) ! ! Write data file of field values within fiducial volume used for fit ! ! The format is the same as the input file. Units are cm and Gauss. ! ! JAC ! 8 March 2016 ! character(*) fname !============================================================ ! Common block of field values from Vector Fields character(80) vfname integer nsize parameter (nsize=1500000) integer nbcalc logical accept(nsize) double precision vfx(nsize) double precision vfbx(nsize) double precision vfebx(nsize) double precision vfy(nsize) double precision vfby(nsize) double precision vfeby(nsize) double precision vfz(nsize) double precision vfbz(nsize) double precision vfebz(nsize) common/vfcom/vfx,vfy,vfz,vfbx,vfby,vfbz, & vfebx,vfeby,vfebz, & nbcalc,accept,vfname !============================================================ character(30) fn double precision zz character(80) string ! Determine step sizes dx=0. dy=0. dz=0. xmax=0. ymax=0. zmax=0. do i=2,nbcalc if(vfx(i-1) /= vfx(i).and.dx == 0.)dx=vfx(i)-vfx(i-1) if(vfy(i-1) /= vfy(i).and.dy == 0.)dy=vfy(i)-vfy(i-1) if(vfz(i-1) /= vfz(i).and.dz == 0.)dz=vfz(i)-vfz(i-1) enddo xmax=0. ymax=0. zmax=0. do i=1,nbcalc if(vfx(i) > xmax)xmax=vfx(i) if(vfy(i) > ymax)ymax=vfy(i) if(vfz(i) > zmax)zmax=vfz(i) enddo nx=1 ny=1 nz=1 if(dx > 0.)nx=int((xmax+dx/2.)/dx)+1 if(dy > 0.)ny=int((ymax+dy/2.)/dy)+1 if(dz > 0.)nz=int((zmax+dz/2.)/dz)+1 write(6,1000)nx,100*dx,ny,100*dy,nz,100*dz 1000 format(' Routine WRITEFV finds',i4,' x bins of width',f4.1,' cm'/ & ' ',i4,' y bins of width',f4.1,' cm'/ & ' ',i4,' z bins of width',f4.1,' cm') z0=0 fn=fname//'_fid.table' print *,' Routine WRITEFV opening file ',fn open(25,file=fn,status='unknown') ! Write number of bins in first line write(25,1100)nx,ny,nz 1100 format(3i10) do ix=1,nx do iy=1,ny do iz=1,nz write(25,2000)100*vfx(i),100*vfy(i),100*vfz(i), & vfbx(i)*10000.,vfbx(i)*10000.,vfbz(i)*10000. 2000 format(6f13.4) enddo enddo enddo close(25) end subroutine readpar(iflag) ! ! Read fit paramter file in BMAD format ! Return nonzero error flag if file not found ! ! JAC ! 1 April 2002 ! ! implicit none !-------------------------------------------------------- ! Fit parameters ! integer NTERMAX parameter (NTERMAX = 600) integer NP parameter (NP = 7*NTERMAX) character(10) pname(NP) integer nterm double precision par(NP),step(NP),pmin(NP),pmax(NP),sigpar(NP) double precision dbx(np),dby(np),dbz(np) common/fitcom/par,step,pmin,pmax,sigpar,pname,nterm common/dbcom/dbx,dby,dbz !-------------------------------------------------------- !----------------------------------------------------------- ! Common block for fit numbering integer numfit common/fitnum/numfit ! !-------------------------------------------------------- ! integer iflag ! character(20) fname integer i,j,kk character(20) abuff integer nt real rpar(7) ! iflag=0 ! ! Get previous fit number call readnum ! write(fname,1000)numfit 1000 format('in/fit',i4.4,'.in') ! open(25,file=fname,err=800,status='old') do i=1,9 read(25,*,err=700,end=700) enddo read(25,2000,err=700,end=700)abuff,nt 2000 format(a13,i10) if(nt <= 0)then write(6,2010)nt 2010 format(' Routine READPAR error: Number of terms in file is',i5) goto 900 elseif(nt > ntermax)then write(6,3000)nt,ntermax 3000 format(' Routine READPAR error: Number of terms in file is',i5/ & ' which exceeds the common block limit of ',i5) goto 900 else do i=1,nt read(25,4000,err=700,end=700)abuff,rpar 4000 format(a10,e13.7,6(1x,e13.7)) ! Coefficients are in Tesla par(1+(i-1)*7)=rpar(1) ! kx,kz are in m^-1 par(2+(i-1)*7)=rpar(2) par(3+(i-1)*7)=rpar(3) ! x0,y0 are in m par(4+(i-1)*7)=rpar(4) par(5+(i-1)*7)=rpar(5) ! Phase is in radians par(6+(i-1)*7)=rpar(6) ! Plane is 0 for X, 1 for Y, 2 for QU and 3 for SQ (see BMAD manual ssect 14.3) par(7+(i-1)*7)=rpar(7) enddo endif close(25) ! ! Write out parameter values write(6,5000)fname,nt 5000 format(' Routine READPAR found file ',a20, & ' Number of terms is',i4) ! write(6,5100)(par(kk),kk=1,nt*7) 5100 format(4f12.6) ! ! Define the number of terms to be used in the fit here nterm=nt ! Now define step,pmin,pmax for C, kx, kz, x0, y0. do i=1,nterm do j=1,3 step(j+(i-1)*7)=abs(par(j+(i-1)*7)) pmin(j+(i-1)*7)=-10*abs(par(j+(i-1)*7)) pmax(j+(i-1)*7)=10*abs(par(j+(i-1)*7)) enddo ! Fix x0,y0 to zero step(4+(i-1)*7)=0. pmin(4+(i-1)*7)=0. pmax(4+(i-1)*7)=0. step(5+(i-1)*7)=0. pmin(5+(i-1)*7)=0. pmax(5+(i-1)*7)=0. ! Fix PHIZ for maximum transverse field at Z=0 step(6+(i-1)*7)=0. pmin(6+(i-1)*7)=0. pmax(6+(i-1)*7)=0. ! Fix family. Do not vary step(7+(i-1)*7)=0. pmin(7+(i-1)*7)=0. pmax(7+(i-1)*7)=3. ! Fix the PHIZ if the field is symmetric around Z=0. Zero if field nonzero at Z=0, otherwise -pi/2. ! step(4+(i-1)*7)=0 enddo ! goto 999 700 continue ! Error during read write(6,7000)fname 7000 format(' Routine READPAR error: Read error or early EOF in file ' & ,a20) goto 900 800 continue ! File not found write(6,8000)fname 8000 format(' Routine READPAR error: File not found -- ',a10) ! 900 continue ! Error condition iflag=1 999 continue ! end subroutine writepar(iflag) ! ! Write fit parameter file in DCS format ! Return nonzero error flag if failure ! ! JAC ! 23 April 2002 ! ! implicit none !============================================================ ! Common block of field values from Vector Fields character(80) vfname integer nsize parameter (nsize=1500000) integer nbcalc logical accept(nsize) double precision vfx(nsize) double precision vfbx(nsize) double precision vfebx(nsize) double precision vfy(nsize) double precision vfby(nsize) double precision vfeby(nsize) double precision vfz(nsize) double precision vfbz(nsize) double precision vfebz(nsize) common/vfcom/vfx,vfy,vfz,vfbx,vfby,vfbz, & vfebx,vfeby,vfebz, & nbcalc,accept,vfname !============================================================ !-------------------------------------------------------- ! Fit parameters ! integer NTERMAX parameter (NTERMAX = 600) integer NP parameter (NP = 7*NTERMAX) character(10) pname(NP) integer nterm double precision par(NP),step(NP),pmin(NP),pmax(NP),sigpar(NP) double precision dbx(np),dby(np),dbz(np) common/fitcom/par,step,pmin,pmax,sigpar,pname,nterm common/dbcom/dbx,dby,dbz !-------------------------------------------------------- !----------------------------------------------------------- ! Common block for fit numbering integer numfit common/fitnum/numfit ! !-------------------------------------------------------- ! integer iflag ! character(80) fname integer i,j,kk character(13) abuff integer nt real rpar(7) ! iflag=0 ! 100 continue numfit=numfit+1 if(numfit > 9999)then print *,' Routine WRITEPAR error condition: NUMFIT=',numfit goto 900 endif print *,' Routine WRITEPAR sets fit number to',numfit ! ! Write final fit parameters to IN directory for use by future job write(fname,1000)numfit 1000 format('in/fit',i4.4,'.in') open(26,file=fname,err=800,status='new') ! call parout(26) ! close(26) write(6,5000)fname 5000 format(' Routine WRITEPAR wrote fit parameters to file ',a20) ! ! Write parameters to file in BMAD format write(fname,1010)numfit 1010 format(i4.4,'.in') fname='bmad/'//vfname(1:index(vfname,' ')-1)//'_'//fname call wbmad(fname) ! ! Write incremented fit number in fitnumber.in call writenum(numfit) ! goto 999 800 continue ! File not found write(6,8000)fname 8000 format(' Routine WRITEPAR error: Could not open file ',a20) goto 100 ! 900 continue ! Error condition iflag=1 999 continue ! end subroutine parout(lunit) ! ! Write fit parameters to LUNIT in DCS format ! implicit none !-------------------------------------------------------- ! Fit parameters ! integer NTERMAX parameter (NTERMAX = 600) integer NP parameter (NP = 7*NTERMAX) character(10) pname(NP) integer nterm double precision par(NP),step(NP),pmin(NP),pmax(NP),sigpar(NP) double precision dbx(np),dby(np),dbz(np) common/fitcom/par,step,pmin,pmax,sigpar,pname,nterm common/dbcom/dbx,dby,dbz !-------------------------------------------------------- ! integer lunit ! character(13) abuff integer i,kk ! abuff=' ' do i=1,9 write(lunit,1010)abuff 1010 format(a13) enddo ! abuff=' n_term = ' write(lunit,1100)abuff,nterm 1100 format(a13,i10) ! do i=1,nterm write(abuff,3000)i 3000 format('term(',i2.2,')=') write(lunit,4000)abuff,(par(7*(i-1)+kk),kk=1,7) ! 4000 format(a13,4f12.6) ! Change to exponential format ! 7/12/2002 4000 format(a10,e13.7,6(1x,e13.7)) enddo ! end subroutine readnum ! ! Read ordinal number of fit from file ! and store it in common block FITNUM ! !----------------------------------------------------------- ! Common block for fit numbering integer numfit common/fitnum/numfit !----------------------------------------------------------- character(20) fname ! fname='fitnumber.in' open(26,file=fname,err=100,status='old') read(26,*)numfit write(6,1000)fname,numfit 1000 format(' Routine READNUM opened file ',a15, & ' Fit number is ',i4) close(26) goto 900 100 continue print *,' Routine READNUM could not open file ',fname numfit=1 ! 900 continue ! end subroutine writenum(numfit) ! ! Write ordinal number of fit from file ! character(20) fname integer numfit ! fname='fitnumber.in' open(26,file=fname,err=100,status='old') write(26,*)numfit close(26) write(6,1000)numfit,fname 1000 format(' Routine WRITENUM incremented fit number to ',i4, & ' in file ',a20) goto 999 100 continue print *,' Routine WRITENUM could not open file ',fname ! 999 continue ! end real function by0x(x) ! use units of gauss and cm for the histograms double precision byfun double precision xx,dzero external byfun data dzero/0./ xx=x/100. by0x=10000.*byfun(xx,dzero,dzero) end real function by0y(y) double precision byfun double precision yy,dzero external byfun data dzero/0./ yy=y/100. by0y=10000.*byfun(dzero,yy,dzero) end real function by0z(z) double precision byfun double precision zz,dzero external byfun data dzero/0./ zz=z/100. by0z=10000.*byfun(dzero,dzero,zz) end subroutine wbmad(fname) ! ! Write fit parameter file in format ! compatible with BMAD tracking code ! ! Return nonzero error flag if failure ! ! JAC ! 1 May 2002 ! ! implicit none !============================================================ ! Common block of field values from Vector Fields character(80) vfname integer nsize parameter (nsize=1500000) integer nbcalc logical accept(nsize) double precision vfx(nsize) double precision vfbx(nsize) double precision vfebx(nsize) double precision vfy(nsize) double precision vfby(nsize) double precision vfeby(nsize) double precision vfz(nsize) double precision vfbz(nsize) double precision vfebz(nsize) common/vfcom/vfx,vfy,vfz,vfbx,vfby,vfbz, & vfebx,vfeby,vfebz, & nbcalc,accept,vfname !============================================================ !-------------------------------------------------------- ! Fit parameters ! integer NTERMAX parameter (NTERMAX = 600) integer NP parameter (NP = 7*NTERMAX) character(10) pname(NP) integer nterm double precision par(NP),step(NP),pmin(NP),pmax(NP),sigpar(NP) double precision dbx(np),dby(np),dbz(np) common/fitcom/par,step,pmin,pmax,sigpar,pname,nterm common/dbcom/dbx,dby,dbz !-------------------------------------------------------- ! real bmadpar(8*NTERMAX) real kx,ky,kz ! character(*) fname ! integer numfit ! character(80) wfname integer i,j,kk character(13) abuff integer nt real rpar(7) ! real zfront real phiz,twopi ! twopi=2*acos(-1.) ! 100 continue ! wfname=fname open(26,file=wfname,err=800,status='new') abuff=' ' do i=1,9 write(26,1010)abuff 1010 format(a13) enddo ! abuff=' n_term = ' write(26,1100)abuff,nterm 1100 format(a13,i10) ! ! Calculate BMADPAR, which is like PAR, except the ky is added do i=1,nterm kx=par(2+7*(i-1)) kz=par(3+7*(i-1)) bmadpar(1+8*(i-1))=par(1+7*(i-1)) bmadpar(2+8*(i-1))=kx bmadpar(4+8*(i-1))=kz bmadpar(5+8*(i-1))=par(4+7*(i-1)) bmadpar(6+8*(i-1))=par(5+7*(i-1)) ! Reference phi to front end of magnet zfront=-vfz(nbcalc) ! print *,' Rtn WBMAD defines front end of magnet at z=',zfront phiz=par(4+7*(i-1))+kz*zfront+twopi bmadpar(7+8*(i-1))=mod(phiz,twopi) if(kx > 0.)then ky=sqrt(kx**2+kz**2) elseif(-kx <= abs(kz))then ky=sqrt(-kx**2+kz**2) elseif(-kx > abs(kz))then ky=sqrt(kx**2-kz**2) endif bmadpar(3+8*(i-1))=ky enddo ! do i=1,nterm-1 write(abuff,3000)i 3000 format('term(',i3.3,')=') write(26,4000)abuff,(bmadpar(8*(i-1)+kk),kk=1,8) ! write(6,4000)abuff,(bmadpar(8*(i-1)+kk),kk=1,8) 4000 format(a10,'{',e13.7,7(',',e13.7),'},&') ! 4000 format(a13,'{',4(f12.6,','),f12.6,'},&') enddo write(abuff,3001)nterm 3001 format('term(',i3.3,')=') write(26,4001)abuff,(bmadpar(8*(nterm-1)+kk),kk=1,8) ! write(6,4001)abuff,(bmadpar(8*(nterm-1)+kk),kk=1,8) 4001 format(a10,'{',e13.7,7(',',e13.7),'}') ! 4001 format(a13,'{',4(f12.6,','),f12.6,'}') close(26) write(6,5000)wfname 5000 format(' Routine WBMAD wrote fit parameters to file ',a40) ! goto 999 800 continue ! File could not be opened write(6,8000)wfname 8000 format(' Routine WBMAD error: Could not open file ',a40) ! 900 continue ! Error condition 999 continue ! end subroutine defuse ! ! Calculate criterion for including a field point in the ! MINUIT chi-squared calculations ! ! This routine assumes that the parameters for calculating ! the fit function have been read in. ! ! 14 August JAC ! implicit none ! !-------------------------------------------------------- double precision bxfun,byfun,bzfun external bxfun external byfun external bzfun !-------------------------------------------------------- !============================================================ ! Common block of field values from Vector Fields character(80) vfname integer nsize parameter (nsize=1500000) integer nbcalc logical accept(nsize) double precision vfx(nsize) double precision vfbx(nsize) double precision vfebx(nsize) double precision vfy(nsize) double precision vfby(nsize) double precision vfeby(nsize) double precision vfz(nsize) double precision vfbz(nsize) double precision vfebz(nsize) common/vfcom/vfx,vfy,vfz,vfbx,vfby,vfbz, & vfebx,vfeby,vfebz, & nbcalc,accept,vfname !============================================================ ! double precision x,y,z,bx,by,bz ! integer i,iuse ! real contcut data contcut/1000./ ! iuse=0 do i=1,nbcalc accept(i)=.false. ! x=vfx(i) y=vfy(i) z=vfz(i) bx=bxfun(x,y,z) by=byfun(x,y,z) bz=bzfun(x,y,z) ! if(abs((vfbx(i)-bx)/vfebx(i)) <= contcut .and. & abs((vfby(i)-by)/vfeby(i)) <= contcut .and. & abs((vfbz(i)-bz)/vfebz(i)) <= contcut) then ! if(abs(vfebx(i)) > 1.e-7 .and. ! . abs(vfeby(i)) > 1.e-7 .and. ! . abs(vfebz(i)) > 1.e-7) then accept(i)=.true. iuse=iuse+1 endif ! enddo ! ! write(6,8100)contcut,iuse,nbcalc ! write(6,8100)iuse,nbcalc 8100 format(' Routine DEFUSE used cut ',f12.1,' to choose ',i6, & ' of ',i6,' field points for the chi2 calculation') ! 8100 format(' Routine DEFUSE used error cut to choose ',i6, ! . ' of ',i6,' field points for the chi2 calculation') ! end subroutine curldiv ! ! Calculate curl and divergence values at each point ! within the fiducial volume ! ! JAC ! 31 October 2002 ! !============================================================ ! Common block of field values from Vector Fields character(80) vfname integer nsize parameter (nsize=1500000) integer nbcalc logical accept(nsize) double precision vfx(nsize) double precision vfbx(nsize) double precision vfebx(nsize) double precision vfy(nsize) double precision vfby(nsize) double precision vfeby(nsize) double precision vfz(nsize) double precision vfbz(nsize) double precision vfebz(nsize) common/vfcom/vfx,vfy,vfz,vfbx,vfby,vfbz, & vfebx,vfeby,vfebz, & nbcalc,accept,vfname !============================================================ !============================================================ ! Common block of curl and divergence values double precision div(nsize) double precision curlx(nsize) double precision curly(nsize) double precision curlz(nsize) common/curlcom/div,curlx,curly,curlz !============================================================ double precision dx,dy,dz double precision x,y,z,x2,y2,z2 double precision bx,by,bz,bx2,by2,bz2 print *, & ' Routine CURLDIV begins calculating divergence and curl values', & ' in units of T/m' ! Determine step sizes and polarity of magnetic field dx=0. dy=0. dz=0. xmax=0. ymax=0. zmax=0. do i=2,nbcalc if(vfx(i-1) /= vfx(i).and.dx == 0.)dx=vfx(i)-vfx(i-1) if(vfy(i-1) /= vfy(i).and.dy == 0.)dy=vfy(i)-vfy(i-1) if(vfz(i-1) /= vfz(i).and.dz == 0.)dz=vfz(i)-vfz(i-1) enddo xmax=0. ymax=0. zmax=0. do i=1,nbcalc if(vfx(i) > xmax)xmax=vfx(i) if(vfy(i) > ymax)ymax=vfy(i) if(vfz(i) > zmax)zmax=vfz(i) enddo nx=int((xmax+dx/2.)/dx)+1 ny=int((ymax+dy/2.)/dy)+1 nz=int((zmax+dz/2.)/dz)+1 do ix=1,nx do iy=1,ny do iz=1,nz i=iz+(iy-1)*nz+(ix-1)*ny*nz ! if(ix == nx)then ii=iz+(iy-1)*nz+(ix-2)*ny*nz div(i)=div(ii) curlx(i)=curlx(ii) curly(i)=curly(ii) curlz(i)=curlz(ii) elseif(iy == ny)then ii=iz+(iy-2)*nz+(ix-1)*ny*nz div(i)=div(ii) curlx(i)=curlx(ii) curly(i)=curly(ii) curlz(i)=curlz(ii) elseif(iz == nz)then ii=iz-1+(iy-1)*nz+(ix-2)*ny*nz div(i)=div(ii) curlx(i)=curlx(ii) curly(i)=curly(ii) curlz(i)=curlz(ii) else ! ix2=iz+(iy-1)*nz+(ix)*ny*nz iy2=iz+(iy)*nz+(ix-1)*ny*nz iz2=iz+1+(iy-1)*nz+(ix-1)*ny*nz ! print *,vfx(ix2),vfx(i) ! print *,vfy(iy2),vfy(i) ! print *,vfz(iz2),vfz(i) ! Divergence div(i)=(vfbx(ix2)-vfbx(i))/dx+ & (vfby(iy2)-vfby(i))/dy+ & (vfbz(iz2)-vfbz(i))/dz ! Curlx curlx(i)=(vfbz(iy2)-vfbz(i))/dy- & (vfby(iz2)-vfby(i))/dz ! Curly curly(i)=(vfbx(iz2)-vfbx(i))/dz- & (vfbz(ix2)-vfbz(i))/dx ! Curlz curlz(i)=(vfby(ix2)-vfby(i))/dx- & (vfbx(iy2)-vfbx(i))/dy endif if(i == 1)write(6,1000) 1000 format(13x,2x,'X(cm)',2x,'Y(cm)',2x,'Z(cm)', & 3x,'BX(T)',4x,'BY(T)',4x,'BZ(T)', & 8x,'Div',6x,'Curlx',6x,'Curly',6x,'Curlz') if(mod(i-1,5000) == 0)then write(6,1100)i,100*vfx(i),100*vfy(i),100*vfz(i), & vfbx(i),vfby(i),vfbz(i), & div(i),curlx(i),curly(i),curlz(i) 1100 format(' Point',i7,3f7.1,3f11.4,4f11.4) endif enddo enddo enddo end subroutine writent ! ! Open file for ntuple, fill ntuple, write out and close file ! ! 7 March 2016 jac ! real hm common/pawc/hm(800000) !------------------------------------------------ integer LUNDAT integer LUNTABLE integer LUNHIS integer LUNTXT common/luncom/LUNDAT,LUNTABLE,LUNHIS,LUNTXT !------------------------------------------------ !----------------------------------------------------------- ! Common block for fit numbering integer numfit common/fitnum/numfit !-------------------------------------------------------- ! character(10) fn ! write(fn,8000)numfit 8000 format('fit',i4.4,'.rz') print *,' Rtn WRITENT opening file ',fn ! Open RZ file for histos and ntuple call hrOPEN(lunhis,'NT',fn,'N',1024,istat) ! print *,' Rtn WRITENT filling ntuple' ! call fillnt ! ! call hldir(' ',' ') call hcdir('//NT',' ') ! call hldir(' ',' ') ! Write histos to file CALL HROUT(0,ICYCLE,' ') ! Close file CALL HREND('NT') close(lunhis) ! print *,' Rtn WRITENT finished writing file ',fn ! 99999 continue end