//
// Daniel Pitzl, Sep 2010
// 12 GeV Zwischenbeschleuniger
//
// madx < zwbl.madx | sed -e "/TWISS/d" > zwbl.log
//
  option, -echo;
  option, -info;

  TITLE, 'HERA-e 27.5 GeV';

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

  Ee   :=  12.0; ! GeV
  epsx := 30e-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 = 0.5e11,
	bunched, kbunch = 2;

  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);
	};
//
//----------------------------------------------------------------------
// magnet definitions
// cp = cqRB
// cp [GeV/c] = 0.3 R[m] B[T]
// take 12 GeV => R = 12/0.3/B
// use N magnets, angle = 360/N
// angle = L/R = 360/N
// => L = R*angle = 12/0.3/B*360/N
//
  B = 1.3;
  rho = Ee/0.3/B;
  N = 248; ! include dispersion suppressors
  ang = 360/N/wt;
  LB = rho*ang;

  print, text = ' ';
  print, text = 'Dipole:';
  value, B;
  value, rho;
  value, N;
  value, ang, ang*wt;
  value, LB;

  BE: RBEND, L = LB, angle = ang;

  QF: QUADRUPOLE, L = 0.70;
  QD: QUADRUPOLE, L = 0.70;

qf->k1 =       0.8658438343 ;
qd->k1 =      -0.8648379176 ;

  SX: SEXTUPOLE,  L = 0.50, k2 =  17.0;
  SY: SEXTUPOLE,  L = 0.50, k2 = -25.0;

! e/p Bx = k2 x y => Bx = p/e k2 x y
! e/p By = 0.5 k2 (x^2-y^2)
! cp/ce = Ee [GeV] / 0.3

  bsx = Ee/0.3 * sx->k2 * 0.03^2; ! 3 cm pole distance
  bsy = Ee/0.3 * sy->k2 * 0.03^2; ! 3 cm pole distance

  print, text = ' ';
  print, text = 'Sextupoles:';
  value, sx->k2;
  value, sy->k2;
  value, bsx, bsy;

  RF: RFCAVITY, L = 0.00, VOLT = 1, FREQ = 499.2, LAG = 0.5; ! for emit

  D040: DRIFT, L =  0.40;
  D030: DRIFT, L =  0.30;
  D020: DRIFT, L =  0.20;
  D010: DRIFT, L =  0.10;
  D005: DRIFT, L =  0.05;
  HK: HKICKER, L = 0.20, KICK = 0.0;
  VK: VKICKER, L = 0.20, KICK = 0.0;

  arcf: line = ( BE, D010, HK, D010, QF, D010, SX, D005, RF, D005
               , BE, D010, VK, D010, QD, D010, SY, D010 );

  arc1: line = ( BE, D010, HK, D010, QF, D010, SX, D005,     D005
               , BE, D010, VK, D010, QD, D010, SY, D010 );
//
//----------------------------------------------------------------------
// just the arc:
//
  use, period = arcf;

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

  match; ! phase advance per fodo arc cell

    vary, name = qf->k1, step = 1E-4;
    vary, name = qd->k1, step = 1E-4;

    global, q1 = 0.20, q2 = 0.20; ! 
    global, q1 = 0.25, q2 = 0.25; ! 
//    global, q1 = 0.30, q2 = 0.30; ! 

    lmdif;

  endmatch;

  print, text = ' ';
  print, text = 'quads:';
  value, qf->k1;
  value, qd->k1;
  dbdx = qf->k1 * Ee/0.3;
  value dbdx;
  bpeak = 0.03*dbdx; ! 3 cm pole distance
  value, bpeak;
  print, text = ' ';

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

  print, text = ' ';
  TWISS, chrom; ! periodic
//
// save beta at start of arc:
//
  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;

  print, text = ' ';
  phasedeg = table(summ,q1)*360;
  value, phasedeg;
  laenge = table(summ,length);
  fdipole = 200*lb / laenge;
  Umfang = n*laenge/2;
  Durchmesser = umfang/pi;
  Radius = durchmesser/2;
  value, laenge, fdipole, Umfang, Durchmesser, Radius;

  print, text = ' ';
  emit; ! 

  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, 'arc cell';
  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;
//
//----------------------------------------------------------------------
// straight section
//
  print, text = '-----------------------------------------------------';
  print, text = 'straight';

  k1qt01 =  0.50;
  k1qt02 = -0.50;

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

  D040: DRIFT, L =  0.40;
  D050: DRIFT, L =  0.50;
  D100: DRIFT, L =  1.00;
  D200: DRIFT, L =  2.00;

  RF06: RFCAVITY, L = 0.60, VOLT = 1.0, FREQ = 499.2, LAG = 0.5;

  rfc: line = (
                d020,
                qt01,
                d020,
                rf06,
                rf06,
                rf06,
                d020,
                qt02,
                d020,
                rf06,
                rf06,
                rf06,
                );

  use, period = rfc;

  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 rfc:
//
  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 rfc
  value, bety0;
  value, alfx0;
  value, alfy0;
  value, dx0;
  value, dy0;
  value, dpx0;
  value, dpy0;

  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;

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

  print, text = ' ';
  TITLE, 'straight from beta';
  PLOT, HAXIS = S, colour=100, interpolate
      , VAXIS = BETX, BETY, vmin = 0, vmax = ymax;
  PLOT, HAXIS = S, colour=100, interpolate, VAXIS = alfX, alfY;

//  stop;
//
//----------------------------------------------------------------------
// from arc into rfc:
//
  print, text = '-----------------------------------------------------';
  print, text = 'arc into rfc';
//
// match magnets:
//
  QM11: QUADRUPOLE, L = 0.7;
  QM12: QUADRUPOLE, L = 0.7;
  QM13: QUADRUPOLE, L = 0.7;
  QM14: QUADRUPOLE, L = 0.7;
//
// need dipoles to suppress dispersion:
//
  BD11: SBEND, L = 1.0;
  BD12: SBEND, L = 1.0;

qm11->k1 =       0.7278417635 ;
qm12->k1 =      -0.7123649977 ;
qm13->k1 =        0.731997593 ;
qm14->k1 =      -0.8274482246 ;
bd11->angle =     0.006881124869 ;
bd12->angle =      0.01644222649 ;

  Dm11: DRIFT, L =  0.40;

  Db11: DRIFT, L =  0.40;

  Dm12: DRIFT, L =  1.40;

  Db12: DRIFT, L =  0.60;

  Dm13: DRIFT, L =  0.90;

  Dm14: DRIFT, L =  1.00;

  Dm15: DRIFT, L =  0.10;

  mtc1: line = (
                dm11,
                qm11,

                db11,
                bd11,

                dm12,
                qm12,

                db12,
                bd12,

                dm13,
                qm13,

                dm14,
                qm14,

                dm15
                );

  m1: marker;

  al1: line = ( 2*arc1, mtc1, m1, 2*rfc );

  use, period = al1;
//
// match arc to rfc
//
  match, beta0 = bet0; ! start in arc
//
    vary, name = qm11->k1, step = 1E-4;
    vary, name = qm12->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;

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

  endmatch;

  print, text = ' ';
  value, qm11->k1;
  value, qm12->k1;
  value, qm13->k1;
  value, qm14->k1;

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

//  stop;

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

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

//  stop;
//
//----------------------------------------------------------------------
// from rfc into arc:
//
  print, text = '-----------------------------------------------------';
  print, text = 'rfc into arc';
//
// match magnets:
//
  QM21: QUADRUPOLE, L = 0.7;
  QM22: QUADRUPOLE, L = 0.7;
  QM23: QUADRUPOLE, L = 0.7;
  QM24: QUADRUPOLE, L = 0.7;
//
// need dipoles to suppress dispersion:
//
  BD21: SBEND, L = 1.0;
  BD22: SBEND, L = 1.0;

qm21->k1 =       0.6715891961 ;
qm22->k1 =      -0.7579156322 ;
qm23->k1 =       0.9479467332 ;
qm24->k1 =      -0.6440567438 ;
bd21->angle =      0.04487910318 ;
bd22->angle =      0.04407262018 ;

  Dm21: DRIFT, L =  0.10;

  Db21: DRIFT, L =  0.40;

  Dm22: DRIFT, L =  0.90;

  Db22: DRIFT, L =  0.90;

  Dm23: DRIFT, L =  1.20;

  Dm24: DRIFT, L =  1.10;

  Dm25: DRIFT, L =  0.10;

  mtc2: line = (
                dm21,
                qm21,

                db21,
                bd21,

                dm22,
                qm22,

                dm23,
                qm23,

                db22,
                bd22,

                dm24,
                qm24,

                dm25
                );

  m2: marker;

  al2: line = ( 2*rfc, mtc2, m2, 2*arc1 );

  use, period = al2;
//
// match rfc into arc
//
  match, betx = betx0, bety = bety0, alfx = alfx0, alfy = alfy0; ! rfc
//
    vary, name = qm21->k1, step = 1E-4;
    vary, name = qm22->k1, step = 1E-4;
    vary, name = qm23->k1, step = 1E-4;
    vary, name = qm24->k1, step = 1E-4;
    vary, name = bd21->angle, step = 1E-4;
    vary, name = bd22->angle, step = 1E-4;

    constraint,
      range = m2, betx = betx5, bety = bety5, alfx = alfx5, alfy = alfy5,dx = dx5, dpx = dpx5;

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

  endmatch;

  print, text = ' ';
  value, qm21->k1;
  value, qm22->k1;
  value, qm23->k1;
  value, qm24->k1;

  value, bd21->angle;
  value, bd22->angle;

//  stop;

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

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

//  stop;
//
//----------------------------------------------------------------------
// ring:
//
  D001: DRIFT, L =  0.01;

  rng0: line = ( 120 * ( arcf, d001 ) );

  use, period = rng0;

  print, text = '-----------------------------------------------------';
  print, text = 'arcs:';

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

  print, text = ' ';
  Umfang = table(summ,length);
  Durchmesser = umfang/pi;
  Radius = durchmesser/2;
  value, Umfang, Durchmesser, Radius;

  print, text = ' ';
  emit; ! 

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

  print, text = '-----------------------------------------------------';
  print, text = 'Ring:';

  rng: line = ( 60*arc1, mtc1, 6*rfc, mtc2
              , 60*arc1, mtc1, 6*rfc, mtc2 );

  use, period = rng;

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

  print, text = ' ';
  Umfang = table(summ,length);
  Durchmesser = umfang/pi;
  Radius = durchmesser/2;
  value, Umfang, Durchmesser, Radius;

  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;

  print, text = ' ';
  TITLE, 'Ring';
  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;

exit;