subroutine design
c---design a graded-beta linac section.  15 parameters are in vv.
c----------------------------------------------------------------------
c
      include 'param_sz.h'
      include 'cfldscom.h'
      include 'constcom.h'
      include 'misccom.h'
      include 'syscom.h'
      include 'ucom.h'
c
      dimension ffc(14)
c--------------------------------------------------------------------------
c*
      if(nn.lt.15)go to 100
c---set convergence criterion to .01 degree
      eps=.001*radian
c---get parameters from vv
      win=vv(1)
      if(win.le.0.)win=w0
      gin=win/erest
      bin=sqrt(gin*(2.+gin))/(1.+gin)
      phin=vv(2)*radian
      e0=vv(3)
      afld=vv(4)
      bfld=vv(5)
      dphi=vv(6)*radian
      nphi=vv(7)
      nmax=vv(8)
      zmax=vv(9)
      wmax=vv(10)
      aper=vv(11)
      nout=vv(12)
      phi0=vv(13)
      nc1=vv(14)
      dwtm=vv(15)
      if(nn.gt.15)ltype=vv(16)
      az=0.
      bz=0.
      cz=0.
      if(nn.gt.16)az=vv(17)
      if(nn.gt.17)bz=vv(18)
      if(nn.gt.18)cz=vv(19)
      attf=0.
      bttf=0.
      cttf=0.
      if(nn.gt.19)attf=vv(20)
      if(nn.gt.20)bttf=vv(21)
      if(nn.gt.21)cttf=vv(22)
      ptotl=0.
      tlen=0.
      if (nc1.gt.0) write(nnout,2)
    2 format(' graded-beta linac section'/
     1 ' cell  length     win    wout    phis      e0       t   phase',
     2 '    ztsq     pcu')
      if (nc1.le.0) write(nnout,3)
    3 format(' graded-beta linac section'/
     1 '  nel  length     win    wout    phis     e0t   phase',
     2 '    ztsq     pcu')
      if(nphi.le.0)dphi=0.
      nc=nc1
      zl=0.
      dw=0.
    5 if(nel.ge.100)go to 120
c---first iteration sequence single cell tank
      pout=phin+dphi+pi
      gw=(win+.5*dw)/erest
      bw=sqrt(gw*(2.+gw))/(1.+gw)*(1.+dphi/pi)
      cl=.5*bw*wavel
c---calculate phase shift to center of cell
   10 ph=phin+.5*pi*bw/bin
c---get transit time factor
      if(nn.le.18)go to 11
      t=attf+bw*(bttf+bw*cttf)
      go to 12
   11 if(ltype.eq.2)call sccfc(bw,ffc)
      if(ltype.eq.7)call dawfc(bw,ffc)
      if(ltype.ne.2.and.ltype.ne.7)call genfc(bw,ffc)
      t=.5*ffc(1)
   12 if(nc.gt.nc1)e0=(afld+bfld*cll/cl)*elast
      dw=.01*e0*t*sin(ph)*cl
      w=win+dw
      g=w/erest
      b=sqrt(g*(2.+g))/(1.+g)
      ph=ph+.5*pi*bw/b
      dp=ph-pout
      if(abs(dp).lt.eps)go to 15
      cl=cl*(1.-dp/(pout-phin))
      bw=2.*cl/wavel
      go to 10
   15 if (nc1.le.0)go to 80
c---continue on using field structure for cell
      dp=pi/72.
   20 b=bin
      w=win
      ph=phin
      dz=cl/36.
      call cellfld(cl,nc)
      if(nc.gt.nc1)e0=(afld+bfld*cll/cl)*elast
      zc=-.5*(cl+dz)
      do 30 i=1,36
      ph=ph+dp*bw/b
      zc=zc+dz
      call cfield(0.,zc,nc,ez,er,bt)
      w=w+.01*e0*ez*sin(ph)*dz
      g=w/erest
      b=sqrt(g*(2.+g))/(1.+g)
   30 ph=ph+dp*bw/b
c---check for convergence
      ep=ph-pout
      if(abs(ep).le.eps)go to 40
      cl=cl/(1.+ep/pi)
      bw=2.*cl/wavel
      go to 20
c---convergence reached.  store parameters
   40 phase=phi0+180.*(nc-nc1)
      phase=amod(phase,360.)
      phis=(ph-pi)/radian
      t=.5*fc(1,nc)
      ztsq=az+bw*(bz+cz*bw)
      pcu=0.
      if(ztsq.gt.0.)pcu=.01*(e0*t)**2*cl/ztsq
      ptotl=ptotl+pcu
      tlen=tlen+cl
      write(nnout,50) nc,cl,win,w,phis,e0,t,phase,ztsq,pcu
   50 format(i5,3f8.3,f8.2,2f8.3,f8.1,2f8.3)
      nel=nel+1
      ntype(nel)=7
      el(4,nel)=phase
      el(5,nel)=.01*e0
      el(6,nel)=nc
      el(7,nel)=dwtm
      phi=phase*radian
      el(9,nel)=sin(phi)
      el(10,nel)=cos(phi)
   60 el(1,nel)=cl
      el(2,nel)=aper
      el(3,nel)=0.
      if(nnout.le.0)go to 70
      if(mod(nc-nc1+1,nout).eq.0)el(3,nel)=1.
   70 elast=e0
      cll=cl
      zloc(nel)=cl
      if(nel.gt.1)zloc(nel)=zloc(nel)+zloc(nel-1)
      nc=nc+1
      nct=nc-nc1
      if(nct.ge.nphi)dphi=0.
      zl=zl+cl
      dw=w-win
      win=w
      bin=b
      phin=ph-pi
      if(nct.ge.nmax)go to 95
      if(zl.ge.zmax)go to 95
      if(w.ge.wmax)go to 95
      go to 5
c---store tank parameters
   80 nel=nel+1
      ntype(nel)=8
      el(4,nel)=.01*e0*t
      el(5,nel)=1.
      el(6,nel)=1.
      phase=phi0+180.*(nc-nc1)
      phase=amod(phase,360.)
      phis=(ph-pi)/radian
      e0t=e0*t
      ztsq=az+bw*(bz+cz*bw)
      pcu=0.
      if(ztsq.gt.0.)pcu=.01*(e0*t)**2*cl/ztsq
      ptotl=ptotl+pcu
      tlen=tlen+cl
      write(nnout,90) nel,cl,win,w,phis,e0t,phase,ztsq,pcu
   90 format(i5,3f8.3,f8.2,f8.3,f8.1,2f8.3)
      el(7,nel)=phase
      el(8,nel)=cl
      el(9,nel)=bw
      el(10,nel)=((wavel/(2.*cl))**2 - 1.)*(twopi/wavel)**2
      go to 60
   95 if(ptotl.gt.0.)write(nnout,96) ptotl
   96 format (' total power =',f8.3,' mw')
      write(nnout,97)tlen
   97 format (' total length of this section is ',f8.2,' cm')
      return
  100 write(ndiag,110) nn
  110 format(' wrong number of parameters on design card.  nn= ',i3)
      call appendparm
      stop ' Abnormal stop design a '
  120 write(ndiag,130)
  130 format('  too many elements being generated in design')
      call appendparm
      stop ' Abnormal stop design b '
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*