subroutine fillbfnt(lun1,lun2,ntid1,ntid2,ntid) * * Read BEAMFIELD data files in LUN1 and LUN2. * These files contain electric field values calculated * by ECLOUD for two different values of beam position. * * Fill header vector BFVEC and ntuple NTID with all data. * * Write 2D (transverse beam-averaging only) and * 3D (include longitudinal beam averaging) averages of fields * and gradients * to vectors which will be written by the macro BEAMFIELD * to files which are formatted * to be read by EFGRAD.KUMAC and GRADTUNES.KUMAC. * * Fill ntuples NTID1 and NTID2 with 2D and 3D averages * * 24 February 2009 JAC * integer lun1,lun2,ntid1,ntid2,ntid * integer i,ii,ib,is integer ix,iy,iz integer nlimit character*80 str real buff1(20),buff2(20) integer ibuff1(20),ibuff2(20) equivalence(ibuff1,buff1) equivalence(ibuff2,buff2) * * Bunch averages real ex1av,ex2av,ey1av,ey2av,dexdxav,deydyav real vec(10) * real zero * data nlimit/50000/ data zero/0./ * * Vector with header info vector bfvec * * 2d and 3d averages to be written to files * for analysis by EFGRAD and GRADTUNES vector bfgrad2d vector bfgrad3d * print *,' FILLBFNT called with units',lun1,lun2, . ' for ntuples ',ntid1,ntid2,ntid *Initialization * pi=acos(-1.) nrec=0 ndata=0 nerr=0 igrid=0 ex1=0. ey1=0. ex2=0. ey2=0. ex1av=0. ey1av=0. ex2av=0. ey2av=0. xysum=0. xyzsum=0. zsum=0. timeav=0. * iz=0 ind2d=0 ind3d=0 * local bunch counter jbunch=0 * * Return here after processing each record 10 continue * Beware: COMIS read gets junk if the record length is too long * Not sure of limit, might be about 100 characters * * if (nrec.eq.0) then read(lun1,9000,err=100,end=900)(buff1(ii),ii=1,12) read(lun2,9000,err=100,end=900)(buff2(ii),ii=1,12) 9000 format(i6,2(e13.6),i6,5(e13.6),2i6,e13.6) nrec=nrec+1 * FILE HEADER: * nbunch,xbeam,ybeam,ngrid,dxgrid,dygrid,xsig,ysig,zsig,nbstep,interspace,bwtnorm * BFVEC: * xbeam1,ybeam1,xbeam2,ybeam2,nbunch,ngrid,dxgrid,dygrid,xsig,ysig,zsig,nbstep,interspace,bwtnorm bfvec(1)=buff1(2) bfvec(2)=buff1(3) bfvec(3)=buff2(2) bfvec(4)=buff2(3) bfvec(5)=ibuff1(1) do ii=1,9 bfvec(ii+5)=buff1(ii+3) enddo bfvec(6)=ibuff1(4) bfvec(12)=ibuff1(10) bfvec(13)=ibuff1(11) x1=bfvec(1) y1=bfvec(2) x2=bfvec(3) y2=bfvec(4) nbunch = int(bfvec(5)) ngrid = int(bfvec(6)) nbstep = int(bfvec(12)) interspace = int(bfvec(13)) xwtsig=2*ngrid/sqrt(12.) ywtsig=2*ngrid/sqrt(12.) zwtsig=(interspace-1)/sqrt(12.) xwtnorm=1./(sqrt(2*pi)*xwtsig) ywtnorm=1./(sqrt(2*pi)*ywtsig) zwtnorm=1./(sqrt(2*pi)*zwtsig) print *,' Step sigmas X/Y/ZWTSIG:', & xwtsig,ywtsig,zwtsig print *,' Normalization factors X/Y/ZWTNORM:', & xwtnorm,ywtnorm,zwtnorm goto 90 elseif (mod(nrec-1,(2*ngrid+1)**2+1).eq.0)then * BUNCH TIME SLICE HEADER: * jb,ib,t,bwt,bweight read(lun1,9010,err=100,end=900)(buff1(ii),ii=1,5) read(lun2,9010,err=100,end=900)(buff2(ii),ii=1,5) 9010 format(2i6,3(e14.6)) nrec=nrec+1 * print *,' Time slice header record.NREC=',nrec ibunch=ibuff1(1) ib=ibuff1(2) time=buff1(3) bwt=buff1(4) * zwt = zwtnorm * exp (-((iz-(interspace-1)/2)**2/(2*zwtsig**2))) zsum=zsum+zwt * c increment time sample counter within bunch iz=iz+1 * if(iz.eq.1)jbunch=jbunch+1 if(ibunch.ne.jbunch)then * New bunch. Reset time slice counter. iz=1 jbunch=ibunch endif * * print *,' Header bunch nr, bunch slice, time, bwt=', * & ibunch,ib,time,bwt * print *,' iz, jbunch=',iz,jbunch * goto 90 else * E-field data record * read(lun1,9020,err=100,end=900)(buff1(ii),ii=1,4) read(lun2,9020,err=100,end=900)(buff2(ii),ii=1,4) 9020 format(4(e15.6)) nrec=nrec+1 * * Exclude records between bunches if (iz.gt.interspace)goto 90 * ndata=ndata+1 igrid=igrid+1 * vec(1)=time*1e-9 vec(2)=ibunch vec(3)=buff1(1) vec(4)=buff1(2) vec(5)=5.7e-12*buff1(3) vec(6)=5.7e-12*buff1(4) vec(7)=buff2(1) vec(8)=buff2(2) vec(9)=5.7e-12*buff2(3) vec(10)=5.7e-12*buff2(4) * if(ndata.gt.nlimit)then write(6,1000)ndata-1 1000 format(' Rtn FILLBFNT reaches fill limit of ',i9, . '. No more ntuple entries will be made.') goto 900 endif * if(ndata.le.nlimit)call hfn(ntid,vec) * * Calculate 2-D transverse and 3-D bunch averages * inner loop is y iy=mod(igrid-1,2*ngrid+1) ix=(igrid-1)/(2*ngrid+1) xwt = xwtnorm * exp (-((ix-ngrid)**2/(2*xwtsig**2))) ywt = ywtnorm * exp (-((iy-ngrid)**2/(2*ywtsig**2))) * * print *,' igrid,ix,iy,iz,xwt,ywt,zwt', * & igrid,ix,iy,iz,xwt,ywt,zwt * ex1=ex1+xwt*ywt*vec(5) ey1=ey1+xwt*ywt*vec(6) ex2=ex2+xwt*ywt*vec(9) ey2=ey2+xwt*ywt*vec(10) xysum=xysum+xwt*ywt * ex1av=ex1av+xwt*ywt*zwt*vec(5) ey1av=ey1av+xwt*ywt*zwt*vec(6) ex2av=ex2av+xwt*ywt*zwt*vec(9) ey2av=ey2av+xwt*ywt*zwt*vec(10) xyzsum=xyzsum+xwt*ywt*zwt * if(igrid.eq.(2*ngrid+1)**2)then * End of time slice * ex1=ex1 * ey1=ey1 * ex2=ex2 * ey2=ey2 dexdx=0. deydy=0. if(x1.ne.x2)dexdx=(ex2-ex1)/(x2-x1) if(y1.ne.y2)deydy=(ey2-ey1)/(y2-y1) exav=(ex1+ex2)/2. eyav=(ey1+ey2)/2. * Correct for normalization sum error if(xysum.gt.0)then ex1=ex1/xysum ey1=ey1/xysum ex2=ex2/xysum ey2=ey2/xysum exav=exav/xysum eyav=eyav/xysum dexdx=dexdx/xysum deydy=deydy/xysum endif * * Store 2D averages ind2d=ind2d+1 bfgrad2d(1,ind2d)=time*1e-9 if(x1.ne.x2)then bfgrad2d(2,ind2d)=ex1 bfgrad2d(4,ind2d)=ex2 bfgrad2d(6,ind2d)=dexdx bfgrad2d(8,ind2d)=exav else bfgrad2d(3,ind2d)=ey1 bfgrad2d(5,ind2d)=ey2 bfgrad2d(7,ind2d)=deydy bfgrad2d(9,ind2d)=eyav endif * write(lun3,5000) * & time*1e-9,ex1,ey1,ex2,ey2, * & dexdx,deydy,exav,eyav * 5000 format(9(1x,e15.5)) * vec(1)=time*1e-9 vec(2)=0 vec(3)=ex1 vec(4)=ey1 vec(5)=ex2 vec(6)=ey2 vec(7)=dexdx vec(8)=deydy vec(9)=exav vec(10)=eyav * print *,' Writing to ntuple at time=',time call hfn(ntid1,vec) * if(jbunch.eq.1.and.iz.eq.1) & print *,' Transverse sum normalization check:',xysum xysum=0 igrid=0 ex1=0. ey1=0. ex2=0. ey2=0. * timeav=timeav+time/interspace * if(iz.eq.interspace)then * End of bunch * Write 3-D averages to mock EFGRAD file * ex1av=ex1av * ey1av=ey1av * ex2av=ex2av * ey2av=ey2av dexdx=0. deydy=0. if(x1.ne.x2)dexdx=(ex2av-ex1av)/(x2-x1) if(y1.ne.y2)deydy=(ey2av-ey1av)/(y2-y1) exav=(ex1av+ex2av)/2. eyav=(ey1av+ey2av)/2. * Correct for normalization sum error if(xyzsum.gt.0)then ex1av=ex1av/xyzsum ey1av=ey1av/xyzsum ex2av=ex2av/xyzsum ey2av=ey2av/xyzsum exav=exav/xyzsum eyav=eyav/xyzsum dexdx=dexdx/xyzsum deydy=deydy/xyzsum endif * Store 3D averages ind3d=ind3d+1 bfgrad3d(1,ind3d)=time*1e-9 if(x1.ne.x2)then bfgrad3d(2,ind3d)=ex1av bfgrad3d(4,ind3d)=ex2av bfgrad3d(6,ind3d)=dexdx bfgrad3d(8,ind3d)=exav else bfgrad3d(3,ind3d)=ey1av bfgrad3d(5,ind3d)=ey2av bfgrad3d(7,ind3d)=deydy bfgrad3d(9,ind3d)=eyav endif * write(lun4,5000) * & time*1e-9,ex1av,ey1av,ex2av,ey2av, * & dexdx,deydy,exav,eyav * * print *,'nrec,iz,jbunch,timeav',nrec,iz,jbunch,timeav vec(1)=timeav*1e-9 vec(2)=0 vec(3)=ex1av vec(4)=ey1av vec(5)=ex2av vec(6)=ey2av vec(7)=dexdx vec(8)=deydy vec(9)=exav vec(10)=eyav call hfn(ntid2,vec) * if(jbunch.eq.1)then print *,' 3D sum normalization check:',xyzsum print *,' Z sum normalization check:',zsum endif * iz=0 xyzsum=0 zsum=0 ex1av=0. ey1av=0. ex2av=0. ey2av=0. timeav=0. endif endif * * * End branch over record type endif * 90 continue * * if(nrec.le.4)then * print *,' Rtn FILLBFNT: Record ',nrec,':',buff1 * endif * * if(nrec.lt.60)goto 10 goto 10 100 continue nrec=nrec+1 nerr=nerr+1 print *, ' Error found in record ',nrec * if(nrec.lt.60)goto 10 goto 10 900 continue write(6,2000)ndata,nrec 2000 format('FILLBFNT finds ',i6, . ' E-field data lines in ',i6,' records') write(6,2010)nerr 2010 format(' ',i6,' lines with format errors found') * end