//
// D. Pitzl, 29.8.2010
// develop new HERA-e ring
// arc cell 4 with 72 deg
// only one IP
// other straight sections with zero dispersion
// SBEND instead RBEND
//
// madx < rng71sb.madx | sed -e "/TWISS/d" > rng71sb.log
//
  option, -echo;
  option, -info;

  TITLE, 'HERA-e 27.5 GeV';

  show, pi; ! MAD knows pi
  wt = 180/pi; ! Winkeltransformation: rad -> deg

  ee   :=  27.5; ! GeV, dE = -108 MeV/turn, eps = 3.4 nm rad
!  ee   :=  45.6; ! GeV, dE = -813 MeV/turn, eps = 9.4 nm rad
  epsx := 3.4e-9; ! pi m rad
  epsy := 0.05*epsx; ! coupling

  beam, particle = positron,
        energy = ee, radiate,
	ex := epsx, ey := epsy,
	sigt := 0.01,
	sige := 1.1e-3,
	npart = 1.5e11,
	bunched, kbunch = 42;

  setplot, lwidth =11;

  coguess, tolerance = 1e-3; ! twiss convergence tolerance, def 1e-6

  macang(ii): macro = {
	ang = table(twiss,angle,ii);
	tlt = table(twiss,tilt,ii);
	vlt = table(twiss,volt,ii);
	};
//
//----------------------------------------------------------------------
// transfer cell
//
  print, text = '-----------------------------------------------------';
  print, text = 'tf2';

  k1qt01 =  0.2055; ! 
  k1qt02 = -0.2045;
/*
table( summ q1 ) =        95.19335073 ;
table( summ q2 ) =        96.71165396 ;
table( summ dq1 ) =         3.28795268 ;
table( summ dq2 ) =        1.535678081 ;
*/

  QT01: QUADRUPOLE, L =  0.5, k1 := k1qt01;
  QT02: QUADRUPOLE, L =  0.5, k1 := k1qt02;

  D00400: DRIFT, L =  0.400;
  D00500: DRIFT, L =  0.500;
  D01000: DRIFT, L =  1.000;
  D02000: DRIFT, L =  2.000;
!  RF12: RFCAVITY, L = 1.20, VOLT = 1.2, FREQ = 499.2, LAG = 0.5;
  RF12: RFCAVITY, L = 1.20, VOLT = 1.2, harmon = 9487, LAG = 0.5;
!  RF12: RFCAVITY, L = 1.20, VOLT = 8, FREQ = 499.2, LAG = 0.5; ! 45 GeV

  tf2: line = (
                d02000,
                qt01,
                d01000,
                d00400,
                rf12,
                d00400,
                d01000,
                qt02,
                d02000
                );

  use, period = tf2;

  print, text = ' ';
  value, qt01->k1;
  value, qt02->k1;

  print, text = ' ';
  TWISS, chrom; ! periodic

  print, text = ' ';
  phasedeg = table(summ,q1)*360;
  value, phasedeg;
//
// save beta at start of tf2:
//
  print, text = ' ';
  print, text = tabstring( twiss, name, 1 );
  s0 = table( twiss, s, 1 );
  betx0 = table( twiss, betx, 1 );
  bety0 = table( twiss, bety, 1 );
  alfx0 = table( twiss, alfx, 1 );
  alfy0 = table( twiss, alfy, 1 );
  dx0 = table( twiss, dx, 1 );
  dy0 = table( twiss, dy, 1 );
  dpx0 = table( twiss, dpx, 1 );
  dpy0 = table( twiss, dpy, 1 );

  value, s0;
  value, betx0; ! start of tf2
  value, bety0;
  value, alfx0;
  value, alfy0;
  value, dx0;
  value, dy0;
  value, dpx0;
  value, dpy0;

  print, text = ' ';
  TITLE, 'tf2 periodic';
  PLOT, HAXIS = S, colour=100, interpolate, VAXIS = BETX, BETY;
//
// tf2 from saved beta:
//
  twiss,betx = betx0, bety = bety0, alfx = alfx0, alfy = alfy0;

  print, text = ' ';
  TITLE, 'tf2 from beta';
  PLOT, HAXIS = S, colour=100, interpolate, VAXIS = BETX, BETY;
  PLOT, HAXIS = S, colour=100, interpolate, VAXIS = alfX, alfY;
//
//----------------------------------------------------------------------
// match focus and tf2
//
  print, text = '-----------------------------------------------------';
  print, text = 'focus and tf2';

  bx0 = 0.03;
  by0 = 0.004;

  value, bx0;
  value, by0;

  TITLE, 'B*x = 3 cm, B*y = 4 mm';
//
// focus magnets:
//
  QF01: QUADRUPOLE, L =  0.8;
  QF02: QUADRUPOLE, L =  0.8;
  QF03: QUADRUPOLE, L =  0.8;
  QF04: QUADRUPOLE, L =  0.4;

qf01->k1 =       -2.191301257 ;
qf02->k1 =        1.552667347 ;
qf03->k1 =      -0.6185143705 ;
qf04->k1 =    4.495804189e-05 ;

  Dfoc: DRIFT, L =  0.500;
  Df01: DRIFT, L =  0.300;
  Df02: DRIFT, L =  0.400;
  Df03: DRIFT, L = 40.000;
  Df04: DRIFT, L = 12.190 - 0.4; ! make optimal to reduce BETy peak

  foc0: line = (
                dfoc,

                qf01,
                df01,

                qf02,
                df02,

                qf03,
                df03,

                qf04,
                df04
                );

  t0: marker;

  ins: line = ( foc0, t0, 8*tf2 );

  use, period = ins;
//
// match focus to transfer
//
  match, betx=bx0, bety=by0, alfx=0, alfy=0;
//
    vary, name = qf01->k1, step = 1E-4;
    vary, name = qf02->k1, step = 1E-4;
    vary, name = qf03->k1, step = 1E-4;
    vary, name = qf04->k1, step = 1E-4;

    constraint,
      range = t0, betx = betx0, bety = bety0, alfx = alfx0, alfy = alfy0;

    lmdif, calls = 5555;
!    migrad;
!    jacobian, calls = 99;

  endmatch;

  print, text = ' ';
  value, qf01->k1;
  value, qf02->k1;
  value, qf03->k1;
  value, qf04->k1;

  print, text = ' ';
  TWISS, chrom, betx=bx0, bety=by0, alfx=0, alfy=0;

  print, text = ' ';
  TITLE, 'B*x = 3 cm, B*y = 4 mm';
  PLOT, HAXIS = S, hmin = 0, hmax = 10, interpolate, colour = 100
      , VAXIS = BETX, BETY;

  print, text = ' ';
  TITLE, 'focus matched tf2';
  PLOT, HAXIS = S, colour = 100, interpolate, VAXIS = BETX, BETY;
  PLOT, HAXIS = S, colour = 100, interpolate, VAXIS = alfX, alfY;
  PLOT, HAXIS = S, colour = 100, interpolate, VAXIS = muX, muY;
  PLOT, HAXIS = S, colour = 100, interpolate
      , VAXIS = dmuX, dmuY, vmin = -50, vmax = 10;

//  stop;
//
//----------------------------------------------------------------------
// old HERA-e BE:
//
  BE: SBEND, L =  9.185, angle =  0.015103813, k1 = 9e-5;

  lbe = BE->L;
  angbe = BE->angle;
  value, lbe, angbe, angbe*wt;
//
//----------------------------------------------------------------------
// double F0D0 arc cell with shorter BE:
//
  angbe04 = 0.5*angbe;!	half the bending for lower emittance
  lbe04 = 0.4*lbe; ! shorter BE, gain space for extra quad
  value, lbe04, angbe04, angbe04*wt;

  angbe04 = 360/365.13*angbe04; ! get 360 degree ring

  BE04    : SBEND, L = lbe04, angle = angbe04;

  QE04    : QUADRUPOLE, L =  0.793, k1 := k1qe04;
  QE05    : QUADRUPOLE, L =  0.793, k1 := k1qe05;

  k1qe04 = -2.64266e-01; ! matched to 2*72 deg
  k1qe05 =  2.64325e-01; ! matched to 2*72 deg

  SL01: SEXTUPOLE,  L =  0.620;
  SK01: SEXTUPOLE,  L =  0.700;

  SL01->k2 =  3.0; ! cures arc dmux chroma for 72 deg
  SK01->k2 = -4.2; ! cures arc dmuy chroma for 72 deg

  SL01->k2 =  5.0; ! cures all dmux chroma for 72 deg
  SK01->k2 = -6.7; ! cures all dmuy chroma for 72 deg

  D00130: DRIFT, L =  0.130;
  D00255: DRIFT, L =  0.255;
  D00082: DRIFT, L =  0.082;
  D00205: DRIFT, L =  0.205;
  D00238: DRIFT, L =  0.238;
  D00212: DRIFT, L =  0.212;
  D00213: DRIFT, L =  0.213;
//
// new arc cell with shorter dipoles (0.4):
// 4 sextupoles
//
  arc_cell4: line = (
              D00130  , ! L =  0.130 at  0.130
              QE04    , ! L =  0.793 at  0.923
              D00255  , ! L =  0.255 at  1.178
              SK01    , ! L =  0.700 at  1.878
              D00213  , ! L =  0.213 at  2.090
              D00205  , ! L =  0.205 at  2.296
              BE04    , ! L =  3.674 at  5.970

              D00238  , ! L =  0.238 at  6.208
              QE05    , ! L =  0.793 at  7.001
              D00255  , ! L =  0.255 at  7.256
              SL01    , ! L =  0.620 at  7.876
              D00205  , ! L =  0.205 at  8.081
              BE04    , ! L =  3.674 at 11.755 = mid cell

              D00130  , ! L =  0.130 at 11.885
              QE04    , ! L =  0.793 at 12.678
              D00255  , ! L =  0.255 at 12.933
              SK01    , ! L =  0.700 at 13.633
              D00213  , ! L =  0.213 at 13.846
              D00205  , ! L =  0.255 at 14.051
              BE04    , ! L =  3.674 at 17.725

              D00238  , ! L =  0.238 at 17.963
              QE05    , ! L =  0.793 at 18.756
              D00255  , ! L =  0.255 at 19.011
              SL01    , ! L =  0.620 at 19.631
              D00205  , ! L =  0.205 at 19.836
              BE04      ! L =  3.674 at 23.510
              );
//
//----------------------------------------------------------------------
// just the new arc cell: double F0D0
//
  use, period = arc_cell4;

  print, text = '-----------------------------------------------------';
  print, text = 'new arc cell';

  match; ! phase advance per double fodo cell: 0.4*360 = 144 degree

    vary,   name = qe04->k1, step = 1E-4;
    vary,   name = qe05->k1, step = 1E-4;

!    global, q1 = 0.33, q2 = 0.33; ! eps arc = 6.0, 60 deg each
    global, q1 = 0.40, q2 = 0.40; ! eps arc = 3.6, 72 deg each
!    global, q1 = 0.60, q2 = 0.60; ! eps arc = 1.2, 108 deg each

    lmdif;

  endmatch;

  savebeta, place = #s, label = bet0; ! start of arc cell, periodic

  print, text = ' ';
  TWISS, chrom;

  print, text = ' ';
  value, beam->ex;

  print, text = ' ';
  value, qe04->k1;
  value, qe05->k1;
  print, text = ' ';

  print, text = ' ';
  phasedeg = 0.5*table(summ,q1)*360; ! per new half-cell
  value, phasedeg;
  print, text = ' ';

  print, text = ' ';
  show, bet0;

  bxmax = table(summ,betxmax);
  bymax = table(summ,betymax);
  if( bxmax > bymax ) {bmax = bxmax;}
  else {bmax = bymax;}
  print, text = ' ';
  value, bxmax, bymax, bmax;
  ymax = bmax+12;
  int imax = ymax/10;
  ymax = 10*imax;
  value, ymax;

  print, text = ' ';
  TITLE, '72 deg arc cell 4';
  PLOT, HAXIS = S, colour = 100, interpolate
      , VAXIS = BETX, BETY, vmin = 0, vmax = ymax;
  PLOT, HAXIS = S, colour = 100, interpolate, VAXIS = DX, DY;
  PLOT, HAXIS = S, colour = 100, interpolate, VAXIS = muX, muY;
  PLOT, HAXIS = S, colour = 100, interpolate, VAXIS = dmuX, dmuY;

//  stop;
//
//----------------------------------------------------------------------
// just the new arc:
//
  arc4: line = ( 4*(arc_cell4) ); ! must be multiple of 4 for periodic

  use, period = arc4;

  print, text = '-----------------------------------------------------';
  print, text = 'new arc4';

  TWISS, chrom; ! periodic

  print, text = ' ';
  value, beam->ex;

  bxmax = table(summ,betxmax);
  bymax = table(summ,betymax);
  if( bxmax > bymax ) {bmax = bxmax;}
  else {bmax = bymax;}
  print, text = ' ';
  value, bxmax, bymax, bmax;
  ymax = bmax+12;
  int imax = ymax/10;
  ymax = 10*imax;
  value, ymax;

  print, text = ' ';
  TITLE, 'new arc4, periodic';
  PLOT, HAXIS = S, colour = 100, interpolate
      , VAXIS = BETX, BETY, vmin = 0, vmax = ymax;
  PLOT, HAXIS = S, colour = 100, interpolate, VAXIS = DX, DY;
//
// save beta at start of arc4:
//
  print, text = ' ';
  print, text = tabstring( twiss, name, 1 );
  s5 = table( twiss, s, 1 );
  betx5 = table( twiss, betx, 1 );
  bety5 = table( twiss, bety, 1 );
  alfx5 = table( twiss, alfx, 1 );
  alfy5 = table( twiss, alfy, 1 );
  dx5 = table( twiss, dx, 1 );
  dy5 = table( twiss, dy, 1 );
  dpx5 = table( twiss, dpx, 1 );
  dpy5 = table( twiss, dpy, 1 );

  value, s5;
  value, betx5;
  value, bety5;
  value, alfx5;
  value, alfy5;
  value, dx5;
  value, dy5;
  value, dpx5;
  value, dpy5;
//
//----------------------------------------------------------------------
// match tf2 to arc4
//
  print, text = '-----------------------------------------------------';
  print, text = 'tf2 and arc4';
//
// match magnets:
//
  QM01: QUADRUPOLE, L = 1.0;
  QM02: QUADRUPOLE, L = 1.0;
  QM031: QUADRUPOLE, L = 0.5;
  QM032: QUADRUPOLE, L = 0.5;
  QM033: QUADRUPOLE, L = 0.5;
  QM04: QUADRUPOLE, L = 0.5;
//
// need dipoles to build up some dispersion:
//
  lng = 2;

  BD01: SBEND, L = lng;
  BD02: SBEND, L = lng;

qm01->k1 =      -0.6846580261 ;
qm02->k1 =       0.4605540664 ;
qm031->k1 =      -0.5879543208 ;
qm032->k1 =      0.04159043452 ;
qm033->k1 =       0.1887099392 ;
qm04->k1 =      0.06438186472 ;

bd01->angle =     0.004672541856 ;
bd02->angle =    -0.009993283822 ;

  Dm01: DRIFT, L =  1.000;
  Dm02: DRIFT, L =  1.500;
  Db01: DRIFT, L =  0.300;
  Dm03: DRIFT, L =  1.000;
  Db02: DRIFT, L =  0.500;
  Dm04: DRIFT, L =  3.000;

  mtc0: line = (
                Qm01,
                Dm01,

                Qm02,
                Dm02,

                Bd01,
                Db01,

                Qm031,
                Dm03,

                Bd01,
                Db01,

                Qm032,
                Dm03,

                Bd01,
                Db01,

                Qm033,
                Dm03,

                Bd02,
                Db02,

                Qm04,
                Dm04,

                Bd02,
                Db02,

                Qm04,
                Dm04
                );

  m0: marker;

  al0: line = ( 4*tf2, mtc0, m0, arc4 );

  use, period = al0;
//
// match tf2 to arc4
//
  match, betx = betx0, bety = bety0, alfx = alfx0, alfy = alfy0; ! tf2
//
    vary, name = qm01->k1, step = 1E-4;
    vary, name = qm02->k1, step = 1E-4;
    vary, name = qm031->k1, step = 1E-4;
    vary, name = qm032->k1, step = 1E-4;
    vary, name = qm033->k1, step = 1E-4;
    vary, name = qm04->k1, step = 1E-4;
    vary, name = bd01->angle, step = 1E-4;
    vary, name = bd02->angle, step = 1E-4;

    constraint,
      range = m0, betx = betx5, bety = bety5, alfx = alfx5, alfy = alfy5,dx = dx5, dpx = dpx5, dy = dy5, dpy = dpy5;

    lmdif, calls = 5555;
!    migrad;
!    jacobian, calls = 99;

  endmatch;

  print, text = ' ';
  value, qm01->k1;
  value, qm02->k1;
  value, qm031->k1;
  value, qm032->k1;
  value, qm033->k1;
  value, qm04->k1;

  value, bd01->angle;
  value, bd02->angle;

//  stop;

  print, text = ' ';
  TWISS, chrom, betx = betx0, bety = bety0, alfx = alfx0, alfy = alfy0;

  print, text = ' ';
  TITLE, 'tf2 mtc arc4';
  PLOT, HAXIS = S, colour=100, interpolate, VAXIS = BETX, BETY;
  PLOT, HAXIS = S, colour=100, interpolate, VAXIS = DX, DY;
  PLOT, HAXIS = S, colour = 100, interpolate, VAXIS = muX, muY;
  PLOT, HAXIS = S, colour = 100, interpolate, VAXIS = dmuX, dmuY;

//  stop;
//
//----------------------------------------------------------------------
// foc to arc:
//
  print, text = '-----------------------------------------------------';
  print, text = 'foc to arc';

  al2: line = ( ins, mtc0, m0, arc4 );

  use, period = al2;

  print, text = ' ';
  TWISS, chrom, betx=bx0, bety=by0, alfx=0, alfy=0;

  print, text = ' ';
  TITLE, 'foc to arc4';
  PLOT, HAXIS = S, colour=100, interpolate, VAXIS = BETX, BETY;
  PLOT, HAXIS = S, colour=100, interpolate, VAXIS = DX, DY;
  PLOT, HAXIS = S, colour = 100, interpolate, VAXIS = muX, muY;
  PLOT, HAXIS = S, colour = 100, interpolate, VAXIS = dmuX, dmuY;

//  stop;
//
//----------------------------------------------------------------------
// match arc4 into tf2:
//
  print, text = '-----------------------------------------------------';
  print, text = 'arc4 and tf2';
//
// match magnets:
//
  QM11: QUADRUPOLE, L = 1.0;
  QM121: QUADRUPOLE, L = 0.5;
  QM122: QUADRUPOLE, L = 0.5;
  QM123: QUADRUPOLE, L = 0.5;
  QM13: QUADRUPOLE, L = 0.5;
  QM14: QUADRUPOLE, L = 0.5;
//
// need dipoles to suppress dispersion:
//
  BD11: SBEND, L = lng;
  BD12: SBEND, L = lng;

qm11->k1 =      -0.2219949656 ;
qm121->k1 =       0.1056412082 ;
qm122->k1 =       0.2705568857 ;
qm123->k1 =       0.0650918943 ;
qm13->k1 =      -0.2420515524 ;
qm14->k1 =        0.245447514 ;
bd11->angle =        0.018206034 ;
bd12->angle =     -0.01909051434 ;

  Dm11: DRIFT, L =  0.500;

  Db11: DRIFT, L =  0.500;
  Dm12: DRIFT, L =  1.000;

  Db12: DRIFT, L =  0.500;
  Dm13: DRIFT, L =  1.000;

  Dm14: DRIFT, L =  1.000;

  mtc1: line = (
                qm11,
                dm11,

                bd11,
                db11,

                qm121,
                dm12,

                bd11,
                db11,

                qm122,
                dm12,

                bd11,
                db11,

                qm123,
                dm12,

                bd12,
                db12,

                qm13,
                dm13,

                bd12,
                db12,

                qm13,
                dm13,

                qm14,
                dm14
                );

  m1: marker;

  al1: line = ( arc4, mtc1, m1, 4*tf2 );

  use, period = al1;
//
// match arc4 to tf2
//
  match, beta0 = bet0; ! start in arc4
//
    vary, name = qm11->k1, step = 1E-4;
    vary, name = qm121->k1, step = 1E-4;
    vary, name = qm122->k1, step = 1E-4;
    vary, name = qm123->k1, step = 1E-4;
    vary, name = qm13->k1, step = 1E-4;
    vary, name = qm14->k1, step = 1E-4;
    vary, name = bd11->angle, step = 1E-4;
    vary, name = bd12->angle, step = 1E-4;

    constraint,
      range = m1, betx = betx0, bety = bety0, alfx = alfx0, alfy = alfy0,dx = 0.0, dpx = 0, dy = 0, dpy = 0;

    lmdif, calls = 5555;
!    migrad;
!    jacobian, calls = 99;

  endmatch;

  print, text = ' ';
  value, qm11->k1;
  value, qm121->k1;
  value, qm122->k1;
  value, qm123->k1;
  value, qm13->k1;
  value, qm14->k1;

  value, bd11->angle;
  value, bd12->angle;

//  stop;

  print, text = ' ';
  TWISS, chrom, beta0 = bet0; ! start in arc4

  print, text = ' ';
  TITLE, 'arc4 into tf2';
  PLOT, HAXIS = S, colour=100, interpolate, VAXIS = BETX, BETY;
  PLOT, HAXIS = S, colour=100, interpolate, VAXIS = DX, DY;

//  stop;
//
//----------------------------------------------------------------------
// quadrant
//
  print, text = '-----------------------------------------------------';
  print, text = 'quadrant';

  qdr: line = ( foc0, 12*tf2, mtc0, 13*arc4, mtc1, 12*tf2 );

  use, period = qdr;

  print, text = ' ';
  TWISS, chrom, betx=bx0, bety=by0, alfx=0, alfy=0; ! from foc0

  setplot, lwidth = 5;

  print, text = ' ';
  TITLE, 'Quadrant';
  PLOT, HAXIS = S, colour=100, VAXIS = BETX, BETY;
  PLOT, HAXIS = S, colour=100, VAXIS = DX, DY;
  PLOT, HAXIS = S, colour = 100, VAXIS = muX, muY;
  PLOT, HAXIS = S, colour = 100, VAXIS = dmuX, dmuY;

//  stop;
//
//----------------------------------------------------------------------
// half
//
  print, text = '-----------------------------------------------------';
  print, text = 'half';

  hlf: line = ( foc0, 12*tf2, mtc0, 13*arc4, mtc1, 16*tf2,
                      16*tf2, mtc0, 13*arc4, mtc1, 12*tf2 );

  use, period = hlf;

  print, text = ' ';
  TWISS, chrom, betx=bx0, bety=by0, alfx=0, alfy=0; ! from foc0

  setplot, lwidth = 5;

  print, text = ' ';
  TITLE, 'half';
  setplot, lwidth = 5;
  PLOT, HAXIS = S, colour=100, VAXIS = BETX, BETY;
  PLOT, HAXIS = S, colour=100, VAXIS = DX, DY;
  PLOT, HAXIS = S, colour = 100, VAXIS = muX, muY;
  PLOT, HAXIS = S, colour = 100, VAXIS = dmuX, dmuY;

//  stop;
//
//----------------------------------------------------------------------
// match tf2 into focus0
//
  print, text = '-----------------------------------------------------';
  print, text = 'tf2 and focus0';
//
// focus magnets:
//
  QF41: QUADRUPOLE, L =  0.4;
  QF42: QUADRUPOLE, L =  0.8;
  QF43: QUADRUPOLE, L =  0.8;
  QF44: QUADRUPOLE, L =  0.8;

qf41->k1 =      0.06126793442 ;
qf42->k1 =      -0.6047164516 ;
qf43->k1 =        1.817738502 ;
qf44->k1 =       -2.533560556 ;

  Df41: DRIFT, L = 12.220 - 0.4; ! make optimal to reduce BETy peak
  Df42: DRIFT, L = 40.000;
  Df43: DRIFT, L =  0.400;
  Df44: DRIFT, L =  0.090;
  Dfc4: DRIFT, L =  0.500;

// 9.98244641e-07

  foc4: line = (
                df41, ! inverted foc0
                qf41,

                df42,
                qf42,

                df43,
                qf43,

                df44,
                qf44,

                dfc4
                );

  f4: marker;

  ip4: line = ( 4*tf2, foc4, f4 );

  use, period = ip4;
//
// match tf2 to focus
//
  match, betx = betx0, bety = bety0, alfx = alfx0, alfy = alfy0; ! tf2
//
    vary, name = qf41->k1, step = 1E-4;
    vary, name = qf42->k1, step = 1E-4;
    vary, name = qf43->k1, step = 1E-4;
    vary, name = qf44->k1, step = 1E-4;

    constraint,
      range = f4, betx=bx0, bety=by0, alfx=0, alfy=0; ! IP

    lmdif, calls = 5555;
!    migrad;
!    jacobian, calls = 99;

  endmatch;

  print, text = ' ';
  value, qf41->k1;
  value, qf42->k1;
  value, qf43->k1;
  value, qf44->k1;

//  stop;

  print, text = ' ';
  TWISS, betx = betx0, bety = bety0, alfx = alfx0, alfy = alfy0; ! tf2

  print, text = ' ';
  TITLE, 'tf2 matched into focus0';
  setplot, lwidth = 11;
  PLOT, HAXIS = S, colour=100, interpolate, VAXIS = BETX, BETY;
  PLOT, HAXIS = S, colour=100, interpolate, VAXIS = dX, dY;

//  stop;
//
//----------------------------------------------------------------------
// dispersion suppressor
//
  print, text = '-----------------------------------------------------';
  print, text = 'dispersion suppressor';

  dsp: line = ( arc4, mtc1, 4*tf2, foc4 );

  use, period = dsp;

  print, text = ' ';
  TWISS, beta0 = bet0; ! arc

  print, text = ' ';
  TITLE, 'dispersion suppressor';
  setplot, lwidth = 11;
  PLOT, HAXIS = S, colour=100, interpolate, VAXIS = BETX, BETY;
  PLOT, HAXIS = S, colour=100, interpolate, VAXIS = dX, dY;

//  stop;
//
//----------------------------------------------------------------------
// ring
//
  print, text = '-----------------------------------------------------';
  print, text = 'ring';

  rng: line = ( 
              foc0, 12*tf2, mtc0, 13*arc4, mtc1, 16*tf2,
                    16*tf2, mtc0, 13*arc4, mtc1, 16*tf2,
                    16*tf2, mtc0, 13*arc4, mtc1, 16*tf2,
                    16*tf2, mtc0, 13*arc4, mtc1, 12*tf2, foc4
              );

  use, period = rng;

  print, text = ' ';
  TWISS, chrom, betx=bx0, bety=by0, alfx=0, alfy=0; ! from foc0

  print, text = ' ';
  emit; ! 

  print, text = ' ';
  value, beam->ex;

  print, text = ' ';
  TITLE, 'Ring';
  setplot, lwidth = 5;
  PLOT, HAXIS = S, hmin = 0, hmax = 6400, colour = 100
      , VAXIS = BETX, BETY;
  PLOT, HAXIS = S, hmin = 0, hmax = 6400, colour = 100
      , VAXIS = DX, DY;
  PLOT, HAXIS = S, hmin = 0, hmax = 6400, colour = 100
      , VAXIS = muX, muY;
  PLOT, HAXIS = S, hmin = 0, hmax = 6400, colour = 100
      , VAXIS = dmuX, dmuY;
  setplot, lwidth = 11;
  PLOT, HAXIS = S, hmin = 6100, hmax = 6350, colour = 100
      , VAXIS = dmuX, dmuY;
//
//----------------------------------------------------------------------
// ring, periodic
//
  print, text = '-----------------------------------------------------';
  print, text = 'ring periodic';

  print, text = ' ';
  TWISS, chrom; ! periodic

  print, text = ' ';
  value, table( twiss, betx, 1 );
  value, table( twiss, bety, 1 );

  tl = table( twiss, tablelength );
  value, tl;
  sumang = 0;
  ii = 0;
  nbend = 0;
  nrf = 0;
  sumvolt = 0;

  while( ii < tl ) {
    ii = ii + 1;
    exec, macang( $ii );
!    value, ang*wt;
    sumang = sumang + ang*cos(tlt);
    sumvolt = sumvolt + vlt;
    if( abs(ang) > 1e-6 ) { nbend = nbend + 1; }
    if( abs(vlt) > 1e-6 ) { nrf = nrf + 1; }
  }
  value, nbend;
  value, sumang*wt;
  diff = sumang-2*pi;
  value, diff;

  print, text = ' ';
  value, nrf;
  value, sumvolt;

  print, text = '/*';
  value, table(summ,q1);
  value, table(summ,q2);
  value, table(summ,dq1);
  value, table(summ,dq2);
  print, text = '*/';

  print, text = ' ';
  emit; ! 

  print, text = ' ';
  value, beam->ex;

//  stop;

  print, text = ' ';
  TITLE, 'Ring, periodic';
  setplot, lwidth = 11;
  PLOT, HAXIS = S, hmin = 0, hmax = 10, interpolate, colour=100
      , VAXIS = BETX, BETY;
  setplot, lwidth = 5;
  PLOT, HAXIS = S, hmin = 0, hmax = 6400, colour = 100
      , VAXIS = BETX, BETY;
  PLOT, HAXIS = S, hmin = 0, hmax = 6400, colour = 100
      , VAXIS = DX, DY;
  PLOT, HAXIS = S, hmin = 0, hmax = 6400, colour = 100
      , VAXIS = muX, muY;
  PLOT, HAXIS = S, hmin = 0, hmax = 6400, colour = 100
      , VAXIS = dmuX, dmuY;
  setplot, lwidth = 11;
  PLOT, HAXIS = S, hmin = 6100, hmax = 6350, colour = 100
      , VAXIS = dmuX, dmuY;

  save, sequence = rng, file = rng71sb.seq8, mad8;

  save, sequence = rng, file = rng71sb.seqx; ! keyword beam crashes...

exit;