subroutine out6dp(iz0,iplace)
c        print emittances
c        iz0 = 0 -- include all good particles
c        iz0 = 1 -- include only good particles with cord(7,i) > 0
c
c        optcon(1)  = 6 -- to invoke OUT6 via the OUTPUT card
c        optcon(2) ne 0 -- use canonical momentum rather than
c                          ordinary momentum in emittance calculation
c        optcon(3) ne 0 -- printer plot of x vs gamma*betax
c                   < 0 -- also print projection plots
c        optcon(4) ne 0 -- plot y vs gamma*betay
c        optcon(5) ne 0 -- plot z vs r
c        optcon(6) ne 0 -- plot z vs gamma*beta
c        optcon(7) ne 0 -- plot z vs del(gamma*beta)/<gamma*beta>
c
c        OUT6 is called with IZ0 = 0 from INPUT if INTYPE = 10,  If want
c        plots at this time, the OUTPUT card must precede the INPUT card
c
c
c-----------------------------------------------------------------------
      implicit double precision (b,g,r,t,x,y,z)
c warning variables en simple precision dans les commons
      real zloc                                        ! systeme
      real gam                                         ! coords
      real beami,brhof,radian,twopi,z0,zlimit          ! const
      save
c
      include 'param_sz.h'
      include 'cfldscom.h'
      include 'constcom.h'
      include 'coordcom.h'
      include 'flagcom.h'
      include 'ncordscom.h'
      include 'outcom.h'
      include 'syscom.h'
      include 'tstepcom.h'
      include 'ucom.h'
      include 'upawcom.h'
      include 'var_char.h'
c
      double precision adp
      common/ddat/ctime,ddat
cbm      common/hist2/xmax,bxmax,ymax,bymax,zminc,zmaxc,rmax,gbmin,
cbm     *             gbmax,gbav,iiz0,ne
cbm      common/hist3/gbxx(imaa),gbyy(imaa)
      dimension gbxx(imaa),gbyy(imaa)
      common/out6/nbphase,nbremesh
      common/outpar/noutp,outpar(0:lmx,10),moutp(10)
c--------------------------------------------------------------------------
c*
c        need double precision for emittance calc
c
      data phold/0./,neold/0/,cp/1./,sp/0./
c        hfac = Compton wavelength times critical magnetic field
      data hfac/1703./
c
      iiz0=iz0
c        copy parameters in case OUT6 not called from OUTPUT
      if(iplace.eq.0) then
         navv=nnout       
         else             
         navv=nav          
       endif              
      do 5 i = 1,lmx
        optcon(i) = outpar(i,7)
    5 continue
      xav = 0.
      xsqav = 0.
      bxav = 0.
      bxsqav = 0.
      xbxav = 0.
      xbxsav = 0.
      yav = 0.
      ysqav = 0.
      byav = 0.
      bysqav = 0.
      ybyav = 0.
      ybysav = 0.
      zav = 0.
      zsqav = 0.
      bzav = 0.
      bzsqav = 0.
      zbzav = 0.
      zbzsav = 0.
      gav = 0.
      gsqav = 0.
      bav = 0.
      bsqav = 0.
      gbav = 0.
      gbsqav = 0.
      zmaxc = 0.
      zminc = 99999.
      adp = 0.
      xmax = 0.
      bxmax = 0.
      ymax = 0.
      bymax = 0.
      rmax = 0.
      gmax = 0.
      gmin = 99999.
      gbmax = 0.
      gbmin = 99999.
      if (ngood.eq.0) go to 110
c######################################################################
      do 100 i = 1,ngood
      ne = cord(7,i)
c        skip if particle is still upstream of first element
      if (ne.eq.0.and.iz0.eq.1) go to 100
      adp = adp + 1.
      x = cord(1,i)
      if (abs(x).gt.xmax) xmax = abs(x)
      xsq = x**2
      y = cord(3,i)
      if (abs(y).gt.ymax) ymax = abs(y)
      ysq = y**2
      r = dsqrt(xsq + ysq)
      if (r.gt.rmax) rmax = r
      z = cord(5,i)
      if (z.gt.zmaxc) zmaxc = z
      if (z.lt.zminc) zminc = z
      ax = 0.
      ay = 0.
      if (optcon(2).eq.0.) go to 20  ! ================================
c        calculate (transverse) vector potential for canonical momentum
c        actually calc eA/mc**2
c           skip if not yet in first element
      if (ne.eq.0) go to 20 ! =========================================
      nt = ntype(ne)
c        skip if not in an rf cell (tho other elements can have nonzero
c                      vector potential...)
      if (nt.ne.7) go to 20
      if (ne.ne.neold) then
        neold = ne
        e0 = el(5,ne)
        s0 = el(9,ne)
        c0 = el(10,ne)
        efac = e0*wavel*(freq/el(4,ne))/twopi/erest
        nc=el(6,ne)
c        spzoff = spzloc(ne) - hcll(nc)
        spzoff = zloc(ne) - hcll(nc)
        if (el(8,ne).gt.0.) spzoff = spzoff + hcll(nc)
        hz = el(11,ne)/hfac/2.
      endif
      ph=wt*radian*el(4,ne)/freq
      if(ph.eq.phold) go to 10
      sp=sin(ph)
      cp=cos(ph)
c        s = sp*c0+cp*s0
        c = cp*c0-sp*s0
      phold=ph
      write(navv,12) e0,s0,c0,wt,c,efac,hz
   12 format (' e0,s0,c0,wt,c,efac,hz =',3f8.4,f8.2,3f8.4)
   10 continue
      if (r.gt.0.) then
        spzc = z - spzoff
        call cfield(r,spzc,nc,ez,er,spbt)
c         rf field
        avec = efac*er*c
        ax = avec*x/r
        ay = avec*y/r
c         static mag field
        ax = ax - y*hz
        ay = ay + x*hz
      endif
   20 continue  ! =====================================================
c        write(navv,11) r,spzc,ez,er,spbt,x,y,ax,cord(2,i),ay,cord(4,i)
   11   format (' r,zc,ez,er,bt =',5f8.4/
     .  ' x,y,ax,bx,ay,by =',6f8.4)
c
c        bx = (canonical momentum)/mc, if 2nd parameter on OUTPUT card ne 0
      bx = cord(2,i)+ ax
      gbxx(i) = bx
      if (abs(bx).gt.bxmax) bxmax = abs(bx)
      bxsq = bx**2
      xbx = x*bx
      xbxsq = xbx**2
      xav = xav + x
      xsqav = xsqav + xsq
      bxav = bxav + bx
      bxsqav = bxsqav + bxsq
      xbxav = xbxav + xbx
      xbxsav = xbxsav + xbxsq
c
      by = cord(4,i) + ay
      gbyy(i) = by
      if (abs(by).gt.bymax) bymax = abs(by)
      bysq = by**2
      yby = y*by
      ybysq = yby**2
      yav = yav + y
      ysqav = ysqav + ysq
      byav = byav + by
      bysqav = bysqav + bysq
      ybyav = ybyav + yby
      ybysav = ybysav + ybysq
c
      zsq = z**2
      bz = cord(6,i)
      bzsq = bz**2
      zbz = z*bz
      zbzsq = zbz**2
      zav = zav + z
      zsqav = zsqav + zsq
      bzav = bzav + bz
      bzsqav = bzsqav + bzsq
      zbzav = zbzav + zbz
      zbzsav = zbzsav + zbzsq
c
c        roundoff trouble extracting beta from gam(i) when gamma near 1
      gb = sqrt(cord(2,i)**2 + cord(4,i)**2 + cord(6,i)**2)
      gbsq = gb**2
      gsq = 1. + gbsq
      g = dsqrt(gsq)
      if (g.gt.gmax) gmax = g
      if (g.lt.gmin) gmin = g
      b = gb/g
c 20/03/2001
      if(b.gt.1.) then
        b = 1.
      endif
      bsq = b**2
      if (gb.gt.gbmax) gbmax = gb
      if (gb.lt.gbmin) gbmin = gb
      gav = gav + g
      gsqav = gsqav + gsq
      bav = bav + b
      bsqav = bsqav + bsq
      gbav = gbav + gb
      gbsqav = gbsqav + gbsq
  100 continue
c#####################################################################
  110 continue
c     1000. rad to mrad
c     1.e6=(1.e3)**2
c --------------------------------- x --------------------------------
      if (adp.gt.0.) then
        xav = xav/adp
        xsqav = xsqav/adp
        bxav = bxav/adp*1000.
        bxsqav = bxsqav/adp*1.e6
        xbxav = xbxav/adp*1000.
        xbxsav = xbxsav/adp*1.e6
      endif
      xrms = 0.
      bxrms = 0.
      xbrms = 0.
      xemit = 0.
      xsqav = xsqav - xav**2
      if (xsqav.gt.0.) xrms = dsqrt(xsqav)
      bxsqav = bxsqav - bxav**2
      if (bxsqav.gt.0.) bxrms = dsqrt(bxsqav)
      temp = xbxsav - xbxav**2
c     if(temp.gt.0) xbrms = sqrt(temp)
      xbxav = xbxav - xav*bxav
      temp = xsqav*bxsqav - xbxav**2
      if (temp.gt.0.) xemit = dsqrt(temp)
c --------------------------------- y --------------------------------
      if (adp.gt.0.) then
        yav = yav/adp
        ysqav = ysqav/adp
        byav = byav/adp*1000.
        bysqav = bysqav/adp*1.e6
        ybyav = ybyav/adp*1000.
        ybysav = ybysav/adp*1.e6
      endif
      yrms = 0.
      byrms = 0.
      ybrms = 0.
      yemit = 0.
      ysqav = ysqav - yav**2
      if (ysqav.gt.0.) yrms = dsqrt(ysqav)
      bysqav = bysqav - byav**2
      if (bysqav.gt.0.) byrms = dsqrt(bysqav)
      temp = ybysav - ybyav**2
c     if(temp.gt.0) ybrms = sqrt(temp)
      ybyav = ybyav - yav*byav
      temp = ysqav*bysqav - ybyav**2
      if (temp.gt.0.) yemit = dsqrt(temp)
c      rrms = 0.5*(xrms + yrms)
c --------------------------------- z --------------------------------
      if (adp.gt.0.) then
        zav = zav/adp
        zsqav = zsqav/adp
        bzav = bzav/adp
        bzsqav = bzsqav/adp
        zbzav = zbzav/adp
        zbzsav = zbzsav/adp
      endif
      zrms = 0.
      zemit = 0.
      zsqav = zsqav - zav**2
      if (zsqav.gt.0.) zrms = dsqrt(zsqav)
      bzsqav = bzsqav - bzav**2
      temp = zbzsav - zbzav**2
      zbzav = zbzav - zav*bzav
      temp = zsqav*bzsqav - zbzav**2
c     10. cm to mm
      if (temp.gt.0.) zemit = 10.*dsqrt(temp)
c --------------------------------- bg --------------------------------
      if (adp.gt.0.) then
        gav = gav/adp
        gsqav = gsqav/adp
        bav = bav/adp
        bsqav = bsqav/adp
        gbav = gbav/adp
        gbsqav = gbsqav/adp
      endif
c      write(navv,199) gav,gsqav,bav,bsqav,gbav,gbsqav
  199 format (' g & b aves ',1pg11.3,5g11.3)
      grms = 0.
      brms = 0.
      gbrms = 0.
      temp = gsqav - gav**2
      if (temp.gt.0.) grms = dsqrt(temp)
      temp = bsqav - bav**2
      if (temp.gt.0.) brms = dsqrt(temp)
      temp = gbsqav - gbav**2
      if (temp.gt.0.) gbrms = dsqrt(temp)
c --------------------------------- write  --------------------------------
      iadp=adp
      write(navv,*)'---------'
      write(navv,190) wt,adp,cord(5,1)  
  190 format (' phase =',f11.2,' degrees, ',f6.0,' active particles',
     1' zrp = ',f10.4,' cm') 
c        convert emittances to mm-mrad (cm-mrad to mm-mrad)
      xemit = xemit*10.
      yemit = yemit*10.
c     zemit=zemit*(360.*1000*erest)/wavel  ! [keV-deg]
c         for now, average x and y quantities
c      xrms = 0.5*(xrms + yrms)
c      bxrms = 0.5*(bxrms + byrms)
c      xbxav = 0.5*(xbxav + ybyav)
c      xemit = 0.5*(xemit + yemit)
      write(navv,200) xrms,bxrms,xbxav,xemit
c changer les f8.4 en f8.2 pour les emit f7.3 en f7.2 pour les gb(x/y)rms
  200 format (' xrms = ',f7.4,' gbxrms = ',f7.2,' xgbxav = ',f8.2,
     . ' xemit = ',f8.2,' mm-mrad')
      xyemit = 0.5*(xemit + yemit)
      write(navv,210) yrms,byrms,ybyav,yemit,xyemit
  210 format (' yrms = ',f7.4,' gbyrms = ',f7.2,' ygbyav = ',f8.2,
     . ' yemit = ',f8.2,' <',f8.2,'>')
      write(navv,220) zminc,zmaxc,zemit
  220 format (' zmin = ',f8.2,' zmax = ',f8.2,' zemit = ',1pg11.3)
      write(navv,230) zav,zrms,gav,grms
  230 format (' zav = ',g12.5,' +- ',g12.5,' gav = ',g12.5,' +- ',g12.5)
      write(navv,231) bav,brms,gbav,gbrms
  231 format (' bav = ',g12.5,' +- ',g12.5,' gbav= ',g12.5,' +- ',g12.5)
      write(nsemit)wt,adp,xrms,bxrms,xbxav,xemit,yrms,byrms,ybyav,yemit,
     *zav,zrms,zminc,zmaxc,zemit
c
c        histogram setup
c
       do 9996 i=1,ngood
       spxx=cord(1,i)
       if ((ne.eq.0).and.(i.eq.1)) then       
          spxxp=0.
          else
          spxxp=cord(2,i)/cord(6,i)
       endif
       spbgx=gbxx(i)
       spyy=cord(3,i)
       if ((ne.eq.0).and.(i.eq.1)) then       
          spyyp=0.
          else
          spyyp=cord(4,i)/cord(6,i)
       endif
       spbgy=gbyy(i)
       spzz=cord(5,i)
       spbgz=cord(6,i)
       phi=wt
       wz=(gam(i)-1.)*erest
c     unit ndes2 at constant phase wt
cbm       write(ndes2)spxx,spxxp,spbgx,spyy,spyyp,spbgy,spzz,
cbm     1 spbgz,phi,wz,ne,i,ngood,i,phizero(i)
       write(ndes2)spxx,spxxp,spbgx,spyy,spyyp,spbgy,spzz,
     1 spbgz,phi,wz,ne,i,iadp,i,phizero(i),ksi1(i),ksi2(i),ksi3(i)
       if(ne.eq.0) then
               spzz=0.
               phi=phizero(i)
c     unit ndes1 end of element z=cte
cbm        write(ndes1)spxx,spxxp,spbgx,spyy,spyyp,spbgy,spzz,
cbm     1  spbgz,phi,wz,ne,i,ngood,i,phizero(i)
       write(ndes1,*)spxx,spxxp,spbgx,spyy,spyyp,spbgy,spzz,
     1 spbgz,phi,wz,ne,i,iadp,i,phizero(i),ksi1(i),ksi2(i),ksi3(i)
        endif
 9996   continue
c        skip if no hists desired
      ihist=optcon(3)+optcon(4)+optcon(5)+optcon(6)+optcon(7)
      if (ihist.eq.0) return      
cbm 18/07/2k      call histogram2
      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*