subroutine srfillnt(lunit,ntid,lunout) * * Read and fill ntuple with synchroton radiation data * * Added length sums and beta averages * 27 Jan 09 jac * integer lunit,ntid,lunout real scutmin,scutmax * real len(10), bx(10), by(10), avlen(10) double precision prate(10), bxprate(10), byprate(10) real pratemax(10), bxprmax(10), byprmax(10) double precision prsig(10), bxprsig(10), byprsig(10) double precision bxpr,bypr * integer nseg(10) character*10 ele(10) data ele/'Dipole','Drift','Wiggler','Quadrupole', . 'Sextupole','Solenoid','Octupole', . 'Non-dipole','Non-drift','Total'/ * integer i character*16 str * character*32 str2 character*2 str2 character*11 str3 character*2 str4 character*100 str5 real buff(9) real rnum real relattr real vec(12) vector rstat vector scut * save sprev,prnormprev * scutmin=scut(1) scutmax=scut(2) * print *,' SRFILLNT called with unit,ntid=',lunit,ntid print *,' and scutmin,scutmax=',scutmin,scutmax * *Initialization * do i=1,6 rstat(i)=0 enddo * totlen=0 do i=1,10 len(i)=0 bx(i)=0 by(i)=0 prate(i)=0 bxprate(i)=0 byprate(i)=0 pratemax(i)=0 bxprmax(i)=0 byprmax(i)=0 prsig(i)=0 bxprsig(i)=0 byprsig(i)=0 enddo * * Assume that the first entry covers * the s region starting at s=0 sprev=0 * nreject=0 nrec=0 ndata=0 10 continue read(lunit,*,err=100,end=900)i,str,buff,str2 * read(lunit,*,err=100,end=900)i,str,buff,str5 * read(lunit,*,err=100,end=900)i,buff,str2,str2,str2,rnum * print *,i,str,buff,str2 * stop * read(lunit,1111,err=100,end=900)i,str,buff,str2 *,buff,str2,str3,str4,relattr *1111 format(i6,a16,f5.3,2e11.3,2e12.4,3f10.3,a27,a5,e11.4) *1111 format(i6,a16,f5.3,2e11.3,e12.4,e12.4,e12.4,3f10.3,a32) * * if(i.le.5)print *,i,str,buff,str2 * * Data record found * ndata=ndata+1 nrec=nrec+1 * s=buff(1) * Use these if spanning s=0 * if(s.gt.768.437/2)then * buff(1)=s-768.437 * else * buff(1)=s * endif * pratio=0 deltas=s-sprev if(deltas.gt.0)then * Calculate photon/sec/m from the photon/sec column plen=buff(3)/deltas totp=buff(5)/deltas * Normalize to 100 mA (6.25e17 e) pr=(buff(6)/deltas)/6.25e17 * * For removing anomalous points if(plen.gt.0.and.prnormprev.gt.0.)then pratio=(pr/plen)/prnormprev endif * sprev=s else totp=0 pr=0 endif * * All these entries were rejected by the PRATIO criterion * if (pr.gt.2e21)then * print *,' prate is too high:',pr * print *,' pratio=',pratio * pr=0 * endif * * print *,' s,deltas,pr:',s,deltas,pr * * Correct "P/len", "Total Power" and "Photons/sec" for the segment length * to obtain Liner Power Density (W/m), Total Power/meter and Photons/m/e buff(3)=plen buff(5)=totp buff(6)=pr * * Define element type indices if(str2(1:2).eq.'SB')then ind=1 elseif(str2(1:2).eq.'WI')then ind=3 elseif(str2(1:2).eq.'QU')then ind=4 elseif(str2(1:2).eq.'SE')then ind=5 elseif(str2(1:2).eq.'SO')then ind=6 elseif(str2(1:2).eq.'OC')then ind=7 else * Drift+RF+Kickers ind=2 endif * do i=1,9 vec(i)=buff(i) enddo vec(10)=ind vec(11)=deltas vec(12)=pratio * * This selection puts an upper limit * on the ratio of the power-normalize rate * to the previous value. The value was * chosen by looking at the four bmad_12wig_20050626a * files * if (pratio.le.20)then * Fill ntuple call hfn(ntid,vec) * Save normalized rate only for accepted entries prnormprev=0 if(plen.gt.0)then prnormprev=pr/plen endif * Store sums if(s.ge.scutmin.and.s.le.scutmax)then totlen=totlen+deltas nseg(ind)=nseg(ind)+1 len(ind)=len(ind)+deltas bx(ind)=bx(ind)+buff(7)*deltas * bx(ind)=bx(ind)+(buff(7)+buff(8))*deltas by(ind)=by(ind)+buff(8)*deltas prate(ind)=prate(ind)+pr*deltas bxpr=pr*buff(7) * bxpr=pr*(buff(7)+buff(8)) bypr=pr*buff(8) bxprate(ind)=bxprate(ind)+bxpr*deltas byprate(ind)=byprate(ind)+bypr*deltas if(pr.gt.pratemax(ind))pratemax(ind)=pr if(bxpr.gt.bxprmax(ind))bxprmax(ind)=bxpr if(bypr.gt.byprmax(ind))byprmax(ind)=bypr prsig(ind)=prsig(ind)+deltas*pr**2 bxprsig(ind)=bxprsig(ind)+deltas*bxpr**2 byprsig(ind)=byprsig(ind)+deltas*bypr**2 if(ind.ne.1)then nseg(8)=nseg(8)+1 len(8)=len(8)+deltas bx(8)=bx(8)+buff(7)*deltas by(8)=by(8)+buff(8)*deltas prate(8)=prate(8)+pr*deltas bxprate(8)=bxprate(8)+bxpr*deltas byprate(8)=byprate(8)+bypr*deltas if(pr.gt.pratemax(8))pratemax(8)=pr if(bxpr.gt.bxprmax(8))bxprmax(8)=bxpr if(bypr.gt.byprmax(8))byprmax(8)=bypr prsig(8)=prsig(8)+deltas*pr**2 bxprsig(8)=bxprsig(8)+deltas*bxpr**2 byprsig(8)=byprsig(8)+deltas*bypr**2 endif if(ind.ne.2)then nseg(9)=nseg(9)+1 len(9)=len(9)+deltas bx(9)=bx(9)+buff(7)*deltas by(9)=by(9)+buff(8)*deltas prate(9)=prate(9)+pr*deltas bxprate(9)=bxprate(9)+bxpr*deltas byprate(9)=byprate(9)+bypr*deltas if(pr.gt.pratemax(9))pratemax(9)=pr if(bxpr.gt.bxprmax(9))bxprmax(9)=bxpr if(bypr.gt.byprmax(9))byprmax(9)=bypr prsig(9)=prsig(9)+deltas*pr**2 bxprsig(9)=bxprsig(9)+deltas*bxpr**2 byprsig(9)=byprsig(9)+deltas*bypr**2 endif nseg(10)=nseg(10)+1 len(10)=len(10)+deltas bx(10)=bx(10)+buff(7)*deltas by(10)=by(10)+buff(8)*deltas prate(10)=prate(10)+pr*deltas bxprate(10)=bxprate(10)+bxpr*deltas byprate(10)=byprate(10)+bypr*deltas if(pr.gt.pratemax(10))pratemax(10)=pr if(bxpr.gt.bxprmax(10))bxprmax(10)=bxpr if(bypr.gt.byprmax(10))byprmax(10)=bypr prsig(10)=prsig(10)+deltas*pr**2 bxprsig(10)=bxprsig(10)+deltas*bxpr**2 byprsig(10)=byprsig(10)+deltas*bypr**2 endif else nreject=nreject+1 endif * if(ndata.eq.1)then * print *,' Rec 1: buff',buff ds=1e7 ssav=s * smin=s smax=s prmin=pr prmax=pr * else * * This means of finding the bin size assumes * adjacent entries will have the smallest difference * if(s.ne.ssav)then if(abs(s-ssav).le.ds)ds=abs(s-ssav) ssav=s ns=ns+1 endif * if(s.lt.smin)smin=s if(s.gt.smax)smax=s if(pr.lt.prmin)prmin=pr if(pr.gt.prmax)prmax=pr * endif * goto 10 100 continue * read(lunit,2000,err=900)abuff 2000 format(1x,a4) nrec=nrec+1 * print *, ' Error found in record ',nrec * keep s value anyway, so deltas will be calculated correctly sprev=buff(1) goto 10 900 continue write(6,1000)ndata,nrec 1000 format('SRFILLNT finds ',i6,' data lines in ',i6,' records') * if(ndata.ge.1) then * is=int(100*(ds+0.005)) ds=is/100. * if(ds.lt.0.0001)ds=0 * if(ns.eq.0)ds=0 * rstat(1)=smin rstat(2)=smax rstat(3)=ds rstat(4)=prmin rstat(5)=prmin print *,' PRATIO rejected ',nreject,' table entries' do i=1,10 if(bx(i).gt.0.)bxprate(i)=(bxprate(i)/bx(i)) if(by(i).gt.0.)byprate(i)=(byprate(i)/by(i)) if(len(i).gt.0)then bx(i)=bx(i)/len(i) by(i)=by(i)/len(i) prate(i)=(prate(i)/len(i)) prsig(i)=(prsig(i)/len(i)) bxprsig(i)=(bxprsig(i)/len(i)) byprsig(i)=(byprsig(i)/len(i)) endif pratemax(i)=pratemax(i) if(bx(i).gt.0.)bxprmax(i)=(bxprmax(i)/bx(i)) if(by(i).gt.0.)byprmax(i)=(byprmax(i)/by(i)) if(bx(i).gt.0.)bxprsig(i)=(bxprsig(i)/(bx(i))**2) if(by(i).gt.0.)byprsig(i)=(byprsig(i)/(by(i))**2) prsig(i)=sqrt(prsig(i)-prate(i)**2) bxprsig(i)=sqrt(bxprsig(i)-bxprate(i)**2) byprsig(i)=sqrt(byprsig(i)-byprate(i)**2) * if(nseg(i).gt.0)then avlen(i)=len(i)/nseg(i) else avlen(i)=0. endif enddo * * Write table of averages write(6,3050) write(lunout,3050) 3050 format('Element',5x,'Nr Seg',3x,'',2x,'Tot Length', & 2x,'Fraction', & 5x,'',5x,'', & 5x,'') do i=1,10 write(6,3100)ele(i),nseg(i),avlen(i),len(i), & 100*len(i)/totlen, & bx(i),by(i), & prate(i) write(lunout,3100)ele(i),nseg(i),avlen(i),len(i), & 100*len(i)/totlen, & bx(i),by(i), & prate(i) 3100 format(a10,2x,i6,2x,f7.3,3x,f8.1,3x,f8.1,'%', & f12.1,f13.1,f15.3) enddo * * Add second table of photon rates, max, avgs, also for beta-weighted rates write(6,3150) write(lunout,3150) 3150 format(/'Element', & 4x,'=',1x,'Max(F)',1x,'RMS(F)', & 1x,'/',1x,'Max(Bx*F/)',1x,'RMS(Bx*F/)', & 1x,'/',1x,'Max(By*F/)',1x,'RMS(By*F/)') do i=1,10 write(6,3200)ele(i), & prate(i),pratemax(i),prsig(i), & bxprate(i),bxprmax(i),bxprsig(i), & byprate(i),byprmax(i),byprsig(i) write(lunout,3200)ele(i), & prate(i),pratemax(i),prsig(i), & bxprate(i),bxprmax(i),bxprsig(i), & byprate(i),byprmax(i),byprsig(i) 3200 format(a10,1x,f9.3,f12.3,f7.3,f10.3,f12.3,f14.3, & f15.3,f12.3,f14.3) enddo * endif * end