Changeset 18 in PSPA for parmelaPSPA

Timestamp:

Sep 24, 2012, 4:31:45 PM (12 years ago)

Author:

lemeur

Message:

corr. bug dans distrir.f

Location:

parmelaPSPA/trunk

Files:

• cread.f (modified) (1 diff)
• distrir.f (modified) (2 diffs)
• input.f (modified) (4 diffs)
• parmelamainv503.f (modified) (1 diff)

parmelaPSPA/trunk/cread.f

@@ -27 +27 @@
      6 'WIGGLER','ALPHAM','STAT','POISSON','SEXTUPOLE'/
       data lfnr,nextr/0,1/  ! lfnr look for next run (or end)
+
       if (jj.eq.100) lfnr=1
    10 continue

parmelaPSPA/trunk/distrir.f

@@ -18 +18 @@
         pas=xmax/nbin
         ntot=0
-c-------
+
+c------- correction Guy Le Meur 24-09-2012
+        do 10 i=1,nbin
+           n(i)=0
+ 10     continue
+
         do 20 i=1,nbin
           if(fac.lt.0.5) then
@@ -33 +38 @@
               n(i)=int(const2*exp(-((i-0.5)*pas)**2/(2*sigsq))
      1             *(i-0.5)*pas*pas+0.5)
-          ntot=ntot+n(i)
+      
+              ntot=ntot+n(i)
               if(ntot.gt.npoints) then
                  n(i)=n(i)-ntot+npoints

parmelaPSPA/trunk/input.f

@@ -100 +100 @@
      * .4987,.4990,.4993,.4995,.4997,.49975,.4998,.49985,.4999,.49995,
      * .5/
+
       do  1 i=1,imaa  ! mcd
       zcath(i)=0.     ! mcd
@@ -598 +599 @@
 c
 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
       if (nn.lt.10) then
         write(nnout,10001)nn
@@ -692 +694 @@
 10007   format (' gun cavity length =',t40,1pg12.3,' cm')
       endif
+
+
 c-----
 c
@@ -707 +711 @@
 c----------------------------- ran = 1 ----
       if(vv(10).eq.1) then
-      do 10008 i = np1,npoints
-c        choose time (phase) of generation, relative to center of the pulse
-10009 continue
-      call rannor(u,v)
-      if (abs(u).gt.tfac) go to 10009
-      phase = u*sigdeg
+         do 10008 i = np1,npoints
+c     choose time (phase) of generation, relative to center of the pulse
+10009       continue
+            call rannor(u,v)
+            if (abs(u).gt.tfac) go to 10009
+            phase = u*sigdeg
 c     choose intial energy at cathode
-10010 continue
-      call rannor(u,v)
-      e = ecath + u*decath
-      if (e.le.0.) go to 10010
-      betcat = sqrt(2.*e/erest)
-      gamma = 1./sqrt(1. - betcat**2)
-c        choose initial momentum, isotropic, but betaz > 0
-      cose = sandom(d)
-      sine = sqrt(1. - cose**2)
-      phi = twopi*sandom(d)
-      bx = betcat*sine*cos(phi)
-      by = betcat*sine*sin(phi)
-      bz = betcat*cose
-      cor(2,i) = gamma*bx
-      cor(4,i) = gamma*by
-      cor(6,i) = gamma*bz
-c        choose (x,y,z)
-10011 continue
-      call rannor(u,v)
-      x0 = u*rsig
-      y0 = v*rsig
-      rsq = x0**2 + y0**2
-      if (rsq.gt.rmaxsq) go to 10011
+10010       continue
+            call rannor(u,v)
+            e = ecath + u*decath
+            if (e.le.0.) go to 10010
+            betcat = sqrt(2.*e/erest)
+            gamma = 1./sqrt(1. - betcat**2)
+c     choose initial momentum, isotropic, but betaz > 0
+            cose = sandom(d)
+            sine = sqrt(1. - cose**2)
+            phi = twopi*sandom(d)
+            bx = betcat*sine*cos(phi)
+            by = betcat*sine*sin(phi)
+            bz = betcat*cose
+            cor(2,i) = gamma*bx
+            cor(4,i) = gamma*by
+            cor(6,i) = gamma*bz
+c     choose (x,y,z)
+10011       continue
+            call rannor(u,v)
+            x0 = u*rsig
+            y0 = v*rsig
+            rsq = x0**2 + y0**2
+            if (rsq.gt.rmaxsq) go to 10011
 c----------------------------------
-      if((spota.eq.0.).or.(spotd.eq.0)) go to 10021
-          anglex10 = atan(x0/spotd)
-          delta = spotd*(1-(cos(spota)/cos(spota+(anglex10))))
-          delta = delta*360./wavel
-          x0= spotd*sin(anglex10)/cos(anglex10+spota)
-          phase=phase+delta
-10021   continue
+            if((spota.eq.0.).or.(spotd.eq.0)) go to 10021
+            anglex10 = atan(x0/spotd)
+            delta = spotd*(1-(cos(spota)/cos(spota+(anglex10))))
+            delta = delta*360./wavel
+            x0= spotd*sin(anglex10)/cos(anglex10+spota)
+            phase=phase+delta
+10021       continue
 c---------------------------------------
-      uph(i)=phase
-      if(uph(i).gt.uphmax)uphmax=uph(i)
-      if(uph(i).lt.uphmin)uphmin=uph(i)      
+            uph(i)=phase
+            if(uph(i).gt.uphmax)uphmax=uph(i)
+            if(uph(i).lt.uphmin)uphmin=uph(i)      
 c---------------------------------------
-      ux0(i)=x0
-      if(ux0(i).gt.ux0max)ux0max=ux0(i)
-      if(ux0(i).lt.ux0min)ux0min=ux0(i)
-      vy0(i)=y0
-      if(vy0(i).gt.vy0max)vy0max=vy0(i)
-      if(vy0(i).lt.vy0min)vy0min=vy0(i)
+            ux0(i)=x0
+            if(ux0(i).gt.ux0max)ux0max=ux0(i)
+            if(ux0(i).lt.ux0min)ux0min=ux0(i)
+            vy0(i)=y0
+            if(vy0(i).gt.vy0max)vy0max=vy0(i)
+            if(vy0(i).lt.vy0min)vy0min=vy0(i)
 c-----
-      z0 = 0.
-      if (ajones.gt.0.) then
-c        calculate the shape of the cathode surface
-        z0 = cathod(rsq)
-        zcath(i) = z0
-      endif
+            z0 = 0.
+            if (ajones.gt.0.) then
+c     calculate the shape of the cathode surface
+               z0 = cathod(rsq)
+               zcath(i) = z0
+            endif
 c-----
-c        shift initial (x,y,z) so that the particle will emerge from
-c        the photocathode at (x0,y0,z0) at time = phase
-      toff = (wtstep + phase)*wavel/360.
-      cor(1,i) = x0 - toff*bx
-      cor(3,i) = y0 - toff*by
-      cor(5,i) = z0 - toff*bz
-10008 continue
-      else  !  vv(10).ne.1
-c----
-      call distriphase(degmax,sigdeg,npoints-1)
-      call distrir(rmax,rsig,npoints-1)
-c----
-      do 10012 i = np1,npoints
-c        choose time (phase) of generation, relative to center of the pulse
-c----
-      u=distp(i-1)
-c----
-      phase=u
-c----
-c        choose intial energy at cathode
-10013 continue
-      call rannor(u,v)
-c----
-      e = ecath + u*decath
-      if (e.le.0.) go to 10013
-      betcat = sqrt(2.*e/erest)
-      gamma = 1./sqrt(1. - betcat**2)
-c        choose initial momentum, isotropic, but betaz > 0
-      cose = sandom(d)
-      sine = sqrt(1. - cose**2)
-      phi = twopi*sandom(d)
-      bx = betcat*sine*cos(phi)
-      by = betcat*sine*sin(phi)
-      bz = betcat*cose
-      cor(2,i) = gamma*bx
-      cor(4,i) = gamma*by
-      cor(6,i) = gamma*bz
-c        choose (x,y,z)
-c----
-c----
-      if(vv(10).eq.2) then 
-10014      continue
-           call rannor(u,v)
-           x0=u*rsig
-           y0=v*rsig
-           rsq=x0**2+y0**2
-           if(rsq.gt.rmaxsq) go to 10014
-      else                              !  vv(10) .ne.2
-           x0 = xr(i-1)
-           y0 = yr(i-1)
-      endif                             !  vv(10) .eq. 2
+c     shift initial (x,y,z) so that the particle will emerge from
+c     the photocathode at (x0,y0,z0) at time = phase
+            toff = (wtstep + phase)*wavel/360.
+            cor(1,i) = x0 - toff*bx
+            cor(3,i) = y0 - toff*by
+            cor(5,i) = z0 - toff*bz
+10008    continue
+      else                      !  vv(10).ne.1
+c---- 
+         call distriphase(degmax,sigdeg,npoints-1)
+         call distrir(rmax,rsig,npoints-1)
+c---- 
+         do 10012 i = np1,npoints
+c     choose time (phase) of generation, relative to center of the pulse
+c---- 
+            u=distp(i-1)
+c---- 
+            phase=u
+c---- 
+c     choose intial energy at cathode
+10013       continue
+            call rannor(u,v)
+c---- 
+            e = ecath + u*decath
+            if (e.le.0.) go to 10013
+            betcat = sqrt(2.*e/erest)
+            gamma = 1./sqrt(1. - betcat**2)
+c     choose initial momentum, isotropic, but betaz > 0
+            cose = sandom(d)
+            sine = sqrt(1. - cose**2)
+            phi = twopi*sandom(d)
+            bx = betcat*sine*cos(phi)
+            by = betcat*sine*sin(phi)
+            bz = betcat*cose
+            cor(2,i) = gamma*bx
+            cor(4,i) = gamma*by
+            cor(6,i) = gamma*bz
+c     choose (x,y,z)
+c---- 
+c---- 
+            if(vv(10).eq.2) then 
+10014          continue
+               call rannor(u,v)
+               x0=u*rsig
+               y0=v*rsig
+               rsq=x0**2+y0**2
+               if(rsq.gt.rmaxsq) go to 10014
+            else                !  vv(10) .ne.2
+               x0 = xr(i-1)
+               y0 = yr(i-1)
+            endif               !  vv(10) .eq. 2
 c----------------------------------
-      if((spota.eq.0.).or.(spotd.eq.0)) go to 10022
-          anglex10 = atan(x0/spotd)
-          delta = spotd*(1-(cos(spota)/cos(spota+(anglex10))))
-          delta = delta*360./wavel
-          x0= spotd*sin(anglex10)/cos(anglex10+spota)
-          phase=phase+delta
-10022   continue
+            if((spota.eq.0.).or.(spotd.eq.0)) go to 10022
+            anglex10 = atan(x0/spotd)
+            delta = spotd*(1-(cos(spota)/cos(spota+(anglex10))))
+            delta = delta*360./wavel
+            x0= spotd*sin(anglex10)/cos(anglex10+spota)
+            phase=phase+delta
+10022       continue
 c---------------------------------------
-      uph(i)=phase
-      if(uph(i).gt.uphmax)uphmax=uph(i)
-      if(uph(i).lt.uphmin)uphmin=uph(i)
-c----
-      ux0(i)=x0
-      if(ux0(i).gt.ux0max)ux0max=ux0(i)
-      if(ux0(i).lt.ux0min)ux0min=ux0(i)
-      vy0(i)=y0
-      if(vy0(i).gt.vy0max)vy0max=vy0(i)
-      if(vy0(i).lt.vy0min)vy0min=vy0(i)
-c----
-      z0 = 0.
-      if (ajones.gt.0.) then
-      rsq = x0**2 + y0**2
-c         calculate the shape of the cathode surface
-          z0 = cathod(rsq)
-          zcath(i) = z0
-      endif
-c        shift initial (x,y,z) so that the particle will emerge from
-c        the photocathode at (x0,y0,z0) at time = phase
-      toff = (wtstep + phase)*wavel/360.
-      cor(1,i) = x0 - toff*bx
-      cor(3,i) = y0 - toff*by
-      cor(5,i) = z0 - toff*bz
-10012 continue
+            uph(i)=phase
+            if(uph(i).gt.uphmax)uphmax=uph(i)
+            if(uph(i).lt.uphmin)uphmin=uph(i)
+c---- 
+            ux0(i)=x0
+            if(ux0(i).gt.ux0max)ux0max=ux0(i)
+            if(ux0(i).lt.ux0min)ux0min=ux0(i)
+            vy0(i)=y0
+            if(vy0(i).gt.vy0max)vy0max=vy0(i)
+            if(vy0(i).lt.vy0min)vy0min=vy0(i)
+c---- 
+            z0 = 0.
+            if (ajones.gt.0.) then
+               rsq = x0**2 + y0**2
+c     calculate the shape of the cathode surface
+               z0 = cathod(rsq)
+               zcath(i) = z0
+            endif
+c     shift initial (x,y,z) so that the particle will emerge from
+c     the photocathode at (x0,y0,z0) at time = phase
+            toff = (wtstep + phase)*wavel/360.
+            cor(1,i) = x0 - toff*bx
+            cor(3,i) = y0 - toff*by
+            cor(5,i) = z0 - toff*bz
+10012    continue
+
 c----------
-      endif ! vv(1).eq.1
+      endif                     ! vv(1).eq.1
 c----------
+
       if(ip.eq.2) then
          open(unit=ninput,file='input10',status='unknown')

parmelaPSPA/trunk/parmelamainv503.f

@@ -346 +346 @@
       open(unit=nin,file='parmin',status='old',
      *     access='sequential',form='formatted')
-      open(unit=nnout,file='parmdebug',status='new',
+      open(unit=nnout,file='parmdebug',status='unknown',
      *        access='sequential',form='formatted')
       call ouvrir