From d5d78b9b8845ff40625de399f650ce85a11ff71e Mon Sep 17 00:00:00 2001 From: zhang <zhang@9a6e40ed-f3a0-4838-9b4a-bf418f78e88d> Date: Mon, 5 Jul 2010 12:14:48 +0000 Subject: [PATCH] ZJ Updqted from tracy2 Tested as working fine --- tracy/tracy/src/max4_lib.cc | 2 - tracy/tracy/src/nsls-ii_lib.cc | 7 +- tracy/tracy/src/physlib.cc | 6301 ++++++++++++++++---------------- tracy/tracy/src/soleillib.cc | 1246 ++++++- tracy/tracy/src/t2elem.cc | 395 +- tracy/tracy/src/t2lat.cc | 32 +- tracy/tracy/src/t2ring.cc | 24 +- tracy/tracy/src/tracy.cc | 1 + 8 files changed, 4708 insertions(+), 3300 deletions(-) diff --git a/tracy/tracy/src/max4_lib.cc b/tracy/tracy/src/max4_lib.cc index cf7f895..cacbb6e 100644 --- a/tracy/tracy/src/max4_lib.cc +++ b/tracy/tracy/src/max4_lib.cc @@ -5,8 +5,6 @@ * * ****************************************************/ - - //copied here form nsls-ii_lib.cc; needed for LoadFieldErr_scl char* get_prm_scl(void) { diff --git a/tracy/tracy/src/nsls-ii_lib.cc b/tracy/tracy/src/nsls-ii_lib.cc index e6bf819..7ba3ac3 100644 --- a/tracy/tracy/src/nsls-ii_lib.cc +++ b/tracy/tracy/src/nsls-ii_lib.cc @@ -3885,7 +3885,7 @@ bool get_nu(const double Ax, const double Ay, const double delta, double x[n_turn], px[n_turn], y[n_turn], py[n_turn]; double nu[2][n_peaks], A[2][n_peaks]; Vector x0; - FILE *fp; + FILE *fp; const bool prt = false; const char file_name[] = "track.out"; @@ -3922,6 +3922,7 @@ bool get_nu(const double Ax, const double Ay, const double delta, } +// tune shift with amplitude void dnu_dA(const double Ax_max, const double Ay_max, const double delta) { bool ok; @@ -3940,7 +3941,7 @@ void dnu_dA(const double Ax_max, const double Ay_max, const double delta) ok = false; fp = file_write("dnu_dAx.out"); - fprintf(fp, "# A_x A_y J_x J_y nu_x nu_y\n"); + fprintf(fp, "# x[mm] y[mm] J_x J_y nu_x nu_y\n"); fprintf(fp, "#\n"); for (i = -n_ampl; i <= n_ampl; i++) { Ax = i*Ax_max/n_ampl; @@ -3963,7 +3964,7 @@ void dnu_dA(const double Ax_max, const double Ay_max, const double delta) ok = false; fp = file_write("dnu_dAy.out"); - fprintf(fp, "# A_x A_y nu_x nu_y\n"); + fprintf(fp, "# x[mm] y[mm] J_x J_y nu_x nu_y\n"); fprintf(fp, "#\n"); for (i = -n_ampl; i <= n_ampl; i++) { Ax = A_min; Ay = i*Ay_max/n_ampl; diff --git a/tracy/tracy/src/physlib.cc b/tracy/tracy/src/physlib.cc index 1b8a9a4..e886c7d 100644 --- a/tracy/tracy/src/physlib.cc +++ b/tracy/tracy/src/physlib.cc @@ -1,3089 +1,3212 @@ -/* Tracy-2 - - J. Bengtsson, CBP, LBL 1990 - 1994 Pascal version - SLS, PSI 1995 - 1997 - M. Boege SLS, PSI 1998 C translation - L. Nadolski SOLEIL 2002 Link to NAFF, Radia field maps - J. Bengtsson NSLS-II, BNL 2004 - - -*/ - - -/**************************/ -/* Routines for printing */ -/**************************/ - -/**** same as asctime in C without the \n at the end****/ -char *asctime2(const struct tm *timeptr) -{ - // terminated with \0. - static char wday_name[7][4] = { - "Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat" - }; - // terminated with \0. - static char mon_name[12][4] = { - "Jan", "Feb", "Mar", "Apr", "May", "Jun", - "Jul", "Aug", "Sep", "Oct", "Nov", "Dec" - }; - static char result[26]; - - sprintf(result, "%.3s %.3s%3d %.2d:%.2d:%.2d %d", - wday_name[timeptr->tm_wday], - mon_name[timeptr->tm_mon], - timeptr->tm_mday, timeptr->tm_hour, - timeptr->tm_min, timeptr->tm_sec, - 1900 + timeptr->tm_year); - return result; -} - -/** Get time and date **/ -struct tm* GetTime() -{ - struct tm *whattime; - /* Get time and date */ - time_t aclock; - time(&aclock); /* Get time in seconds */ - whattime = localtime(&aclock); /* Convert time to struct */ - return whattime; -} - - -void printglob(void) -{ - printf("\n***************************************************************" - "***************\n"); - printf("\n"); - printf(" dPcommon = %9.3e dPparticle = %9.3e" - " Energy [GeV] = %.3f\n", - globval.dPcommon, globval.dPparticle, globval.Energy); - printf(" MaxAmplx [m] = %9.3e MaxAmply [m] = %9.3e" - " RFAccept [%%] = \xB1%4.2f\n", - Cell[0].maxampl[X_][0], Cell[0].maxampl[Y_][0], - globval.delta_RF*1e2); - printf(" MatMeth = %s ", globval.MatMeth ? "TRUE " : "FALSE"); - printf(" Cavity_On = %s ", globval.Cavity_on ? "TRUE " : "FALSE"); - printf(" Radiation_On = %s \n", - globval.radiation ? "TRUE " : "FALSE"); - printf(" bpm = %3d qt = %3d ", - globval.bpm, globval.qt); - printf(" Chambre_On = %s \n", status.chambre ? "TRUE " : "FALSE"); - printf(" hcorr = %3d vcorr = %3d\n\n", - globval.hcorr, globval.vcorr); - printf(" alphac = %22.16e\n", globval.Alphac); - printf(" nux = %19.16f nuz = %19.16f", - globval.TotalTune[X_], globval.TotalTune[Y_]); - if (globval.Cavity_on) - printf(" omega = %11.9f\n", globval.Omega); - else { - printf("\n"); - printf(" ksix = %10.6f ksiz = %10.6f\n", - globval.Chrom[X_], globval.Chrom[Y_]); - } - printf("\n"); - printf(" OneTurn matrix:\n"); - printf("\n"); - prtmat(ss_dim, globval.OneTurnMat); - fflush(stdout); -} - - -double int_curly_H1(long int n) -{ - /* Integration with Simpson's Rule */ - - double curly_H; - Vector2 alpha[3], beta[3], nu[3], eta[3], etap[3]; - - // only works for matrix style calculations - get_twiss3(n, alpha, beta, nu, eta, etap); - - curly_H = (get_curly_H(alpha[0][X_], beta[0][X_], eta[0][X_], etap[0][X_]) - +4.0*get_curly_H(alpha[1][X_], beta[1][X_], - eta[1][X_], etap[1][X_]) - + get_curly_H(alpha[2][X_], beta[2][X_], eta[2][X_], etap[2][X_])) - /6.0; - - return curly_H; -} - - -void prt_sigma(void) -{ - long int i; - double code = 0.0; - FILE *outf; - - outf = file_write("../out/sigma.out"); - - fprintf(outf, "# name s sqrt(sx) sqrt(sx') sqrt(sy) sqrt(sy')\n"); - fprintf(outf, "# [m] [mm] [mrad] [mm] [mrad]\n"); - fprintf(outf, "#\n"); - - for (i = 0; i <= globval.Cell_nLoc; i++) { - switch (Cell[i].Elem.Pkind) { - case drift: - code = 0.0; - break; - case Mpole: - if (Cell[i].Elem.M->Pirho != 0) - code = 0.5; - else if (Cell[i].Elem.M->PBpar[Quad+HOMmax] != 0) - code = sgn(Cell[i].Elem.M->PBpar[Quad+HOMmax]); - else if (Cell[i].Elem.M->PBpar[Sext+HOMmax] != 0) - code = 1.5*sgn(Cell[i].Elem.M->PBpar[Sext+HOMmax]); - else if (Cell[i].Fnum == globval.bpm) - code = 2.0; - else - code = 0.0; - break; - default: - code = 0.0; - break; - } - fprintf(outf, "%4ld %.*s %6.2f %4.1f %9.3e %9.3e %9.3e %9.3e\n", - i, SymbolLength, Cell[i].Elem.PName, Cell[i].S, code, - 1e3*sqrt(Cell[i].sigma[x_][x_]), - 1e3*sqrt(fabs(Cell[i].sigma[x_][px_])), - 1e3*sqrt(Cell[i].sigma[y_][y_]), - 1e3*sqrt(fabs(Cell[i].sigma[y_][py_]))); - } - - fclose(outf); -} - - -void recalc_S(void) -{ - long int k; - double S_tot; - - S_tot = 0.0; - for (k = 0; k <= globval.Cell_nLoc; k++) { - S_tot += Cell[k].Elem.PL; Cell[k].S = S_tot; - } -} - - -bool getcod(double dP, long &lastpos) -{ - return GetCOD(globval.CODimax, globval.CODeps, dP, lastpos); -} - - -void get_twiss3(long int loc, - Vector2 alpha[], Vector2 beta[], Vector2 nu[], - Vector2 eta[], Vector2 etap[]) -{ - /* Get Twiss functions at magnet entrance-, center-, and exit. */ - - long int j, k; - Vector2 dnu; - Matrix Ascr0, Ascr1; - - // Lattice functions at the magnet entrance - for (k = 0; k <= 1; k++) { - alpha[0][k] = Cell[loc-1].Alpha[k]; beta[0][k] = Cell[loc-1].Beta[k]; - nu[0][k] = Cell[loc-1].Nu[k]; - eta[0][k] = Cell[loc-1].Eta[k]; etap[0][k] = Cell[loc-1].Etap[k]; - } - - UnitMat(5, Ascr0); - for (k = 0; k <= 1; k++) { - Ascr0[2*k][2*k] = sqrt(Cell[loc-1].Beta[k]); Ascr0[2*k][2*k+1] = 0.0; - Ascr0[2*k+1][2*k] = -Cell[loc-1].Alpha[k]/Ascr0[2*k][2*k]; - Ascr0[2*k+1][2*k+1] = 1/Ascr0[2*k][2*k]; - } - Ascr0[0][4] = eta[0][X_]; Ascr0[1][4] = etap[0][X_]; - Ascr0[2][4] = eta[1][Y_]; Ascr0[3][4] = etap[1][Y_]; - CopyMat(5, Ascr0, Ascr1); MulLMat(5, Cell[loc].Elem.M->AU55, Ascr1); - dnu[0] = (atan(Ascr1[0][1]/Ascr1[0][0])-atan(Ascr0[0][1]/Ascr0[0][0])) - /(2.0*M_PI); - dnu[1] = (atan(Ascr1[2][3]/Ascr1[2][2])-atan(Ascr0[2][3]/Ascr0[2][2])) - /(2.0*M_PI); - - // Lattice functions at the magnet center - for (k = 0; k <= 1; k++) { - alpha[1][k] = -Ascr1[2*k][2*k]*Ascr1[2*k+1][2*k] - - Ascr1[2*k][2*k+1]*Ascr1[2*k+1][2*k+1]; - beta[1][k] = pow(Ascr1[2*k][2*k], 2.0) + pow(Ascr1[2*k][2*k+1], 2); - nu[1][k] = nu[0][k] + dnu[k]; - j = 2*k + 1; eta[1][k] = Ascr1[j-1][4]; etap[1][k] = Ascr1[j][4]; - } - - CopyMat(5, Ascr1, Ascr0); MulLMat(5, Cell[loc].Elem.M->AD55, Ascr1); - dnu[0] = (atan(Ascr1[0][1]/Ascr1[0][0])-atan(Ascr0[0][1]/Ascr0[0][0])) - /(2.0*M_PI); - dnu[1] = (atan(Ascr1[2][3]/Ascr1[2][2])-atan(Ascr0[2][3]/Ascr0[2][2])) - /(2.0*M_PI); - - // Lattice functions at the magnet exit - for (k = 0; k <= 1; k++) { - alpha[2][k] = -Ascr1[2*k][2*k]*Ascr1[2*k+1][2*k] - - Ascr1[2*k][2*k+1]*Ascr1[2*k+1][2*k+1]; - beta[2][k] = pow(Ascr1[2*k][2*k], 2.0) + pow(Ascr1[2*k][2*k+1], 2); - nu[2][k] = nu[1][k] + dnu[k]; - j = 2*k + 1; eta[2][k] = Ascr1[j-1][4]; etap[2][k] = Ascr1[j][4]; - } -} - - -void getabn(Vector2 &alpha, Vector2 &beta, Vector2 &nu) -{ - Vector2 gamma; - Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu); -} - - -void TraceABN(long i0, long i1, const Vector2 &alpha, const Vector2 &beta, - const Vector2 &eta, const Vector2 &etap, const double dP) -{ - long i, j; - double sb; - ss_vect<tps> Ascr; - - UnitMat(6, globval.Ascr); - for (i = 1; i <= 2; i++) { - sb = sqrt(beta[i-1]); j = i*2 - 1; - globval.Ascr[j-1][j-1] = sb; globval.Ascr[j-1][j] = 0.0; - globval.Ascr[j][j - 1] = -(alpha[i-1]/sb); globval.Ascr[j][j] = 1/sb; - } - globval.Ascr[0][4] = eta[0]; globval.Ascr[1][4] = etap[0]; - globval.Ascr[2][4] = eta[1]; globval.Ascr[3][4] = etap[1]; - - for (i = 0; i < 6; i++) - globval.CODvect[i] = 0.0; - globval.CODvect[4] = dP; - - if (globval.MatMeth) - Cell_Twiss_M(i0, i1, globval.Ascr, false, false, dP); - else { - for (i = 0; i <= 5; i++) { - Ascr[i] = tps(globval.CODvect[i]); - for (j = 0; j <= 5; j++) - Ascr[i] += globval.Ascr[i][j]*tps(0.0, j+1); - } - Cell_Twiss(i0, i1, Ascr, false, false, dP); - } - -} - - -void FitTune(long qf, long qd, double nux, double nuy) -{ - long i; - iVector2 nq = {0,0}; - Vector2 nu = {0.0, 0.0}; - fitvect qfbuf, qdbuf; - - /* Get elements for the first quadrupole family */ - nq[X_] = GetnKid(qf); - for (i = 1; i <= nq[X_]; i++) - qfbuf[i-1] = Elem_GetPos(qf, i); - - /* Get elements for the second quadrupole family */ - nq[Y_] = GetnKid(qd); - for (i = 1; i <= nq[Y_]; i++) - qdbuf[i - 1] = Elem_GetPos(qd, i); - - nu[X_] = nux; nu[Y_] = nuy; - - /* fit tunes */ - Ring_Fittune(nu, nueps, nq, qfbuf, qdbuf, nudkL, nuimax); -} - - -void FitChrom(long sf, long sd, double ksix, double ksiy) -{ - long i; - iVector2 ns = {0,0}; - fitvect sfbuf, sdbuf; - Vector2 ksi ={0.0, 0.0}; - - /* Get elements for the first sextupole family */ - ns[X_] = GetnKid(sf); - for (i = 1; i <= ns[X_]; i++) - sfbuf[i-1] = Elem_GetPos(sf, i); - - /* Get elements for the second sextupole family */ - ns[Y_] = GetnKid(sd); - for (i = 1; i <= ns[Y_]; i++) - sdbuf[i-1] = Elem_GetPos(sd, i); - - ksi[X_] = ksix; ksi[Y_] = ksiy; - - /* Fit chromaticities */ - /* Ring_Fitchrom(ksi, ksieps, ns, sfbuf, sdbuf, 1.0, 1);*/ - Ring_Fitchrom(ksi, ksieps, ns, sfbuf, sdbuf, ksidkpL, ksiimax); -} - - -void FitDisp(long q, long pos, double eta) -{ - long i, nq; - fitvect qbuf; - - /* Get elements for the quadrupole family */ - nq = GetnKid(q); - for (i = 1; i <= nq; i++) - qbuf[i-1] = Elem_GetPos(q, i); - - Ring_FitDisp(pos, eta, dispeps, nq, qbuf, dispdkL, dispimax); -} - - -#define nfloq 4 - -void getfloqs(Vector &x) -{ - // Transform to Floquet space - LinTrans(nfloq, globval.Ascrinv, x); -} - -#undef nfloq - - -#define ntrack 4 - -// 4D tracking in normal or Floquet space over nmax turns - -void track(const char *file_name, - double ic1, double ic2, double ic3, double ic4, double dp, - long int nmax, long int &lastn, long int &lastpos, int floqs, - double f_rf) -{ - /* Single particle tracking around closed orbit: - - Output floqs - - Phase Space 0 [x, px, y, py, delta, ct] - Floquet Space 1 [x^, px^, y^, py^, delta, ct] - Action-Angle Variables 2 [2Jx, phx, 2Jy, phiy, delta, ct] - - */ - long int i; - double twoJx, twoJy, phix, phiy, scl_1 = 1.0, scl_2 = 1.0; - Vector x0, x1, x2, xf; - FILE *outf; - - bool prt = false; - - if ((floqs == 0)) { - scl_1 = 1e3; scl_2 = 1e3; - x0[x_] = ic1; x0[px_] = ic2; x0[y_] = ic3; x0[py_] = ic4; - } else if ((floqs == 1)) { - scl_1 = 1.0; scl_2 = 1.0; - x0[x_] = ic1; x0[px_] = ic2; x0[y_] = ic3; x0[py_] = ic4; - LinTrans(4, globval.Ascr, x0); - } else if (floqs == 2) { - scl_1 = 1e6; scl_2 = 1.0; - x0[x_] = sqrt(ic1)*cos(ic2); x0[px_] = -sqrt(ic1)*sin(ic2); - x0[y_] = sqrt(ic3)*cos(ic4); x0[py_] = -sqrt(ic3)*sin(ic4); - LinTrans(4, globval.Ascr, x0); - } - - outf = file_write(file_name); - fprintf(outf, "# Tracking with TRACY"); - getcod(dp, lastpos); - if (floqs == 0) - fprintf(outf, "\n"); - else if (floqs == 1) { - Ring_GetTwiss(false, dp); - fprintf(outf, "# (Floquet space)\n"); - } else if (floqs == 2) { - Ring_GetTwiss(false, dp); - fprintf(outf, "# (Action-Angle variables)\n"); - } - fprintf(outf, "#\n"); - fprintf(outf, "#%3d%6ld% .1E% .1E% .1E% .1E% 7.5f% 7.5f\n", - 1, nmax, 1e0, 1e0, 0e0, 0e0, globval.TotalTune[0], - globval.TotalTune[1]); - if (floqs == 0) { - fprintf(outf, "# N x p_x y p_y"); - fprintf(outf, " delta cdt\n"); - fprintf(outf, "# [mm] [mrad]" - " [mm] [mrad]"); - } else if (floqs == 1) { - fprintf(outf, "# N x^ px^ y^ py^"); - fprintf(outf, " delta cdt"); - fprintf(outf, "# " - " "); - } else if (floqs == 2) { - fprintf(outf, "# N 2Jx phi_x 2Jy phi_y"); - fprintf(outf, " delta cdt\n"); - fprintf(outf, "# " - " "); - } - if (f_rf == 0.0){ - fprintf(outf, " [%%] [mm]\n"); - fprintf(outf, "#\n"); - fprintf(outf, "%4d %23.16e %23.16e %23.16e %23.16e %23.16e %23.16e\n", - 0, scl_1*ic1, scl_2*ic2, scl_1*ic3, scl_2*ic4, 1e2*dp, - 1e3*globval.CODvect[ct_]); - } else { - fprintf(outf, " [%%] [deg]\n"); - fprintf(outf, "#\n"); - fprintf(outf, "%4d %23.16e %23.16e %23.16e %23.16e %23.16e %23.16e\n", - 0, scl_1*ic1, scl_2*ic2, scl_1*ic3, scl_2*ic4, 1e2*dp, - 2.0*f_rf*180.0*globval.CODvect[ct_]/c0); - } - x2[x_] = x0[x_] + globval.CODvect[x_]; - x2[px_] = x0[px_] + globval.CODvect[px_]; - x2[y_] = x0[y_] + globval.CODvect[y_]; - x2[py_] = x0[py_] + globval.CODvect[py_]; - if (globval.Cavity_on) { - x2[delta_] = dp + globval.CODvect[delta_]; x2[ct_] = globval.CODvect[ct_]; - } else { - x2[delta_] = dp; x2[ct_] = 0.0; - } - - lastn = 0; - - if (prt) { - printf("\n"); - printf("track:\n"); - printf("%4ld %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", - lastn, 1e3*x2[x_], 1e3*x2[px_], 1e3*x2[y_], 1e3*x2[py_], - 1e2*x2[delta_], 1e3*x2[ct_]); - } - - if (globval.MatMeth) Cell_Concat(dp); - - do { - (lastn)++; - for (i = 0; i < nv_; i++) - x1[i] = x2[i]; - - if (globval.MatMeth) { - Cell_fPass(x2, lastpos); - } else { - Cell_Pass(0, globval.Cell_nLoc, x2, lastpos); - } - - for (i = x_; i <= py_; i++) - xf[i] = x2[i] - globval.CODvect[i]; - - for (i = delta_; i <= ct_; i++) - if (globval.Cavity_on && (i != ct_)) - xf[i] = x2[i] - globval.CODvect[i]; - else - xf[i] = x2[i]; - - if (floqs == 1) - getfloqs(xf); - else if (floqs == 2) { - getfloqs(xf); - twoJx = pow(xf[x_], 2) + pow(xf[px_], 2); - twoJy = pow(xf[y_], 2) + pow(xf[py_], 2); - phix = atan2(xf[px_], xf[x_]); - phiy = atan2(xf[py_], xf[y_]); - xf[x_] = twoJx; xf[px_] = phix; xf[y_] = twoJy; xf[py_] = phiy; - } - if (f_rf == 0.0) - fprintf(outf, - "%4ld %23.16le %23.16le %23.16le %23.16le %23.16le %23.16le\n", - lastn, scl_1*xf[0], scl_2*xf[1], scl_1*xf[2], scl_2*xf[3], - 1e2*xf[4], 1e3*xf[5]); - else - fprintf(outf, - "%4ld %23.16le %23.16le %23.16le %23.16le %23.16le %23.16le\n", - lastn, scl_1*xf[0], scl_2*xf[1], scl_1*xf[2], scl_2*xf[3], - 1e2*xf[4], 2.0*f_rf*180.0*xf[5]/c0); - } while ((lastn != nmax) && (lastpos == globval.Cell_nLoc)); - - if (globval.MatMeth) - Cell_Pass(0, globval.Cell_nLoc, x1, lastpos); - - fclose(outf); -} - -#undef ntrack - - -#define step 0.1 -#define px 0.0 -#define py 0.0 -void track_(double r, struct LOC_getdynap *LINK) -{ - long i, lastn, lastpos; - Vector x; - - x[0] = r * cos(LINK->phi); - x[1] = px; - x[2] = r * sin(LINK->phi); - x[3] = py; - x[4] = LINK->delta; - x[5] = 0.0; - /* transform to phase space */ - if (LINK->floqs) { - LinTrans(5, globval.Ascr, x); - } - for (i = 0; i <= 3; i++) - x[i] += globval.CODvect[i]; - lastn = 0; - do { - lastn++; - if (globval.MatMeth) { - Cell_fPass(x, lastpos); - } else { - Cell_Pass(0, globval.Cell_nLoc, x, lastpos); - } - } while (lastn != LINK->nturn && lastpos == globval.Cell_nLoc); - LINK->lost = (lastn != LINK->nturn); -} -#undef step -#undef px -#undef py - - -/****************************************************************************/ -/* void Trac(double x, double px, double y, double py, double dp, double ctau, - long nmax, long pos, long &lastn, long &lastpos, FILE *outf1) - - Purpose: - Single particle tracking w/ respect to the chamber centrum - - Input: - x, px, y, py 4 transverses coordinates - ctau c*tau - dp energy offset - nmax number of turns - pos starting position for tracking - aperture global physical aperture - - Output: - lastn last n (should be nmax if not lost) - lastpos last position in the ring - - Return: - none - - Global variables: - globval - - specific functions: - Cell_Pass - - Comments: - Absolute TRACKING with respect to the center of the vacuum vessel - BUG: last printout is wrong because not at pos but at the end of - the ring - 26/04/03 print output for phase space is for position pos now - 01/12/03 tracking from 0 to pos -1L instead of 0 to pos - (wrong observation point) - -****************************************************************************/ -void Trac(double x, double px, double y, double py, double dp, double ctau, - long nmax, long pos, long &lastn, long &lastpos, FILE *outf1) -{ - bool lostF; /* Lost particle Flag */ - Vector x1; /* tracking coordinates */ - Vector2 aperture; - - /* Compute closed orbit : usefull if insertion devices */ - - aperture[0] = 1e0; - aperture[1] = 1e0; - - x1[0] = x; x1[1] = px; - x1[2] = y; x1[3] = py; - x1[4] =dp; x1[5] = ctau; - - lastn = 0L; - lostF = true; - - (lastpos)=pos; - if(trace) fprintf(outf1, "\n"); - fprintf(outf1, "%6ld %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e \n", lastn, - x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]); - Cell_Pass(pos -1L, globval.Cell_nLoc, x1, lastpos); - - if (lastpos == globval.Cell_nLoc) - do { - (lastn)++; - Cell_Pass(0L,pos-1L, x1, lastpos); - if(!trace) { - fprintf(outf1, "%6ld %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e \n", - lastn, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]); - } - if (lastpos == pos-1L) Cell_Pass(pos-1L,globval.Cell_nLoc, x1, lastpos); - } - while (((lastn) < nmax) && ((lastpos) == globval.Cell_nLoc)); - - if (lastpos == globval.Cell_nLoc) Cell_Pass(0L,pos, x1, lastpos); - - if (lastpos != pos) { - printf("Trac: Particle lost \n"); - fprintf(stdout, "turn:%6ld plane: %1d" - " %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", - lastn, status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]); - } -} - -/****************************************************************************/ -/*bool chk_if_lost(double x0, double y0, double delta, - long int nturn, bool floqs) - - Purpose: - Binary search for dynamical aperture in Floquet space. - - Input: - none - - Output: - none - - Return: - none - - Global variables: - px_0, py_0 - - Specific functions: - chk_if_lost - - Comments: - none - -****************************************************************************/ - -#define nfloq 4 -bool chk_if_lost(double x0, double y0, double delta, - long int nturn, bool floqs) -{ - long int i, lastn, lastpos; - Vector x; - - bool prt = false; - - x[x_] = x0; x[px_] = px_0; x[y_] = y0; x[py_] = py_0; - x[delta_] = delta; x[ct_] = 0.0; - if (floqs) - // transform to phase space - LinTrans(nfloq, globval.Ascr, x); - for (i = 0; i <= 3; i++) - x[i] += globval.CODvect[i]; - - lastn = 0; - if (prt) { - printf("\n"); - printf("chk_if_lost:\n"); - printf("%4ld %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", - lastn, 1e3*x[x_], 1e3*x[px_], 1e3*x[y_], 1e3*x[py_], - 1e2*x[delta_], 1e3*x[ct_]); - } - do { - lastn++; - if (globval.MatMeth) - Cell_fPass(x, lastpos); - else - Cell_Pass(0, globval.Cell_nLoc, x, lastpos); - if (prt) - printf("%4ld %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", - lastn, 1e3*x[x_], 1e3*x[px_], 1e3*x[y_], 1e3*x[py_], - 1e2*x[delta_], 1e3*x[ct_]); - } while ((lastn != nturn) && (lastpos == globval.Cell_nLoc)); - return(lastn != nturn); -} -#undef nfloq - -/****************************************************************************/ -/* void getdynap(double *r, double phi, double delta, double eps, - int nturn, bool floqs) - - Purpose: - Binary search for dynamical aperture in Floquet space. - - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - chk_if_lost - - Comments: - none - -****************************************************************************/ -void getdynap(double &r, double phi, double delta, double eps, - int nturn, bool floqs) -{ - /* Determine dynamical aperture by binary search. */ - double rmin = 0.0, rmax = r; - - const bool prt = false; - const double r_reset = 1e-3, r0 = 10e-3; - - if (prt) printf("\n"); - - if (globval.MatMeth) Cell_Concat(delta); - while (!chk_if_lost(rmax*cos(phi), rmax*sin(phi), delta, nturn, floqs)) { - if (rmax < r_reset) rmax = r0; - rmax *= 2.0; - } - while (rmax-rmin >= eps) { - r = rmin + (rmax-rmin)/2.0; - if (prt) printf("getdynap: %6.3f %6.3f %6.3f\n", - 1e3*rmin, 1e3*rmax, 1e3*r); - if (! chk_if_lost(r*cos(phi), r*sin(phi), delta, nturn, floqs) ) - rmin = r; - else - rmax = r; - if (rmin > rmax) { - printf("getdynap: rmin > rmax\n"); - exit_(0); - } - - } - r = rmin; -} - - - -/****************************************************************************/ -/* void getcsAscr(void) - - Purpose: - Get Courant-Snyder Ascr - - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - none - -****************************************************************************/ -void getcsAscr(void) -{ - long i, j; - double phi; - Matrix R; - - UnitMat(6, R); - for (i = 1; i <= 2; i++) { - phi = -atan2(globval.Ascr[i * 2 - 2][i * 2 - 1], globval.Ascr[i * 2 - 2] - [i * 2 - 2]); - R[i * 2 - 2][i * 2 - 2] = cos(phi); - R[i * 2 - 1][i * 2 - 1] = R[i * 2 - 2][i * 2 - 2]; - R[i * 2 - 2][i * 2 - 1] = sin(phi); - R[i * 2 - 1][i * 2 - 2] = -R[i * 2 - 2][i * 2 - 1]; - } - MulRMat(6, globval.Ascr, R); - for (i = 1; i <= 2; i++) { - if (globval.Ascr[i * 2 - 2][i * 2 - 2] < 0.0) { - for (j = 0; j <= 5; j++) { - globval.Ascr[j][i * 2 - 2] = -globval.Ascr[j][i * 2 - 2]; - globval.Ascr[j][i * 2 - 1] = -globval.Ascr[j][i * 2 - 1]; - } - } - } - if (!InvMat(6, globval.Ascrinv)) - printf(" *** Ascr is singular\n"); -} - - -/****************************************************************************/ -/* void dynap(double r, double delta, double eps, int npoint, int nturn, - double x[], double y[], bool floqs, bool print) - - Purpose: - Determine the dynamical aperture by tracking using polar coordinates, - and sampling in phase. - Assumes mid-plane symmetry - - Input: - r initial guess - delta off momentum energy - eps precision for binary search - npoint sample number for phase coordinate - nturn number of turn for computing da - floqs true means Floquet space - print true means Print out to the screen - - Output: - x[] horizontal dynamics aperture - y[] vertical dynamics aperture - - Return: - none - - Global variables: - none - - Specific functions: - getdynap - - Comments: - none - -****************************************************************************/ -void dynap(FILE *fp, double r, const double delta, - const double eps, const int npoint, const int nturn, - double x[], double y[], const bool floqs, const bool print) - -{ - /* Determine the dynamical aperture by tracking. - Assumes mid-plane symmetry. */ - - long int i, lastpos; - double phi, x_min, x_max, y_min, y_max; - - getcod(delta, lastpos); - if (floqs) { - Ring_GetTwiss(false, delta); - if (print) { - printf("\n"); - printf("Dynamical Aperture (Floquet space):\n"); - printf(" x^ y^\n"); - printf("\n"); - } - fprintf(fp, "# Dynamical Aperture (Floquet space):\n"); - fprintf(fp, "# x^ y^\n"); - fprintf(fp, "#\n"); - } else { - if (print) { - printf("\n"); - printf("Dynamical Aperture:\n"); - printf(" x y\n"); - printf(" [mm] [mm]\n"); - printf("\n"); - } - fprintf(fp, "# Dynamical Aperture:\n"); - fprintf(fp, "# x y\n"); - fprintf(fp, "# [mm] [mm]\n"); - fprintf(fp, "#\n"); - } - - x_min = 0.0; x_max = 0.0; y_min = 0.0; y_max = 0.0; - for (i = 0; i < npoint; i++) { - phi = i*M_PI/(npoint-1); - if (i == 0) - phi = 1e-3; - else if (i == npoint-1) - phi -= 1e-3; - getdynap(r, phi, delta, eps, nturn, floqs); - x[i] = r*cos(phi); y[i] = r*sin(phi); - x_min = min(x[i], x_min); x_max = max(x[i], x_max); - y_min = min(y[i], y_min); y_max = max(y[i], y_max); - if (!floqs) { - if (print) - printf(" %6.2f %6.2f\n", 1e3*x[i], 1e3*y[i]); - fprintf(fp, " %6.2f %6.2f\n", 1e3*x[i], 1e3*y[i]); - } else { - if (print) - printf(" %10.3e %10.3e\n", x[i], y[i]); - fprintf(fp, " %10.3e %10.3e\n", x[i], y[i]); - } - fflush(fp); - } - - if (print) { - printf("\n"); - printf(" x^: %6.2f - %5.2f y^: %6.2f - %5.2f mm\n", - 1e3*x_min, 1e3*x_max, 1e3*y_min, 1e3*y_max); - } -} - -/****************************************************************************/ -/* double get_aper(int n, double x[], double y[]) - - Purpose: - - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - none - -****************************************************************************/ -double get_aper(int n, double x[], double y[]) -{ - int i; - double A; - - A = 0.0; - for (i = 2; i <= n; i++) - A += x[i-2]*y[i-1] - x[i-1]*y[i-2]; - A += x[n-1]*y[0] - x[0]*y[n-1]; -// x2 from mid-plane symmetry - A = fabs(A); -// printf("\n"); -// printf(" Dyn. Aper.: %5.1f mm^2\n", 1e6*A); - return(A); -} - - -void GetTrack(const char *file_name, - long *n, double x[], double px[], double y[], double py[]) -{ - int k; - char line[200]; - FILE *inf; - - inf = file_read(file_name); - - do { - fgets(line, 200, inf); - } while (strstr(line, "#") != NULL); - - // skip initial conditions - fgets(line, 200, inf); - - do { - sscanf(line, "%d", &k); - sscanf(line, "%*d %lf %lf %lf %lf", &x[k-1], &px[k-1], &y[k-1], &py[k-1]); - } while (fgets(line, 200, inf) != NULL); - - *n = k; - - fclose(inf); -} - - -/****************************************************************************/ -/* void Getj(long n, double *x, double *px, double *y, double *py) - - Purpose: - Calculates the linear invariant - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - none - -****************************************************************************/ -void Getj(long n, double *x, double *px, double *y, double *py) -{ - long int i; - - for (i = 0; i < n; i++) { - x[i] = (pow(x[i], 2)+pow(px[i], 2))/2.0; - y[i] = (pow(y[i], 2)+pow(py[i], 2))/2.0; - } -} - -/****************************************************************************/ -/* double GetArg(double x, double px, double nu) - - Purpose: - get argument of x - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - 17/07/03 use M_PI instead of pi - -****************************************************************************/ -double GetArg(double x, double px, double nu) -{ - double phi, val; - - phi = GetAngle(x, px); - if (phi < 0.0) - phi += 2.0 * M_PI; - val = phi + Fract(nu) * 2.0 * M_PI; - if (val < 0.0) - val += 2.0 * M_PI; - return val; -} - -/****************************************************************************/ -/* void GetPhi(long n, double *x, double *px, double *y, double *py) - - Purpose: - get linear phases - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - none - -****************************************************************************/ -void GetPhi(long n, double *x, double *px, double *y, double *py) -{ - /* Calculates the linear phase */ - long i; - - for (i = 1; i <= n; i++) { - x[i - 1] = GetArg(x[i - 1], px[i - 1], i * globval.TotalTune[0]); - y[i - 1] = GetArg(y[i - 1], py[i - 1], i * globval.TotalTune[1]); - } -} - - -/*********************************/ -/* Routines for Fourier analysis */ -/*********************************/ - -void Sinfft(int n, double xr[]) -{ - /* DFT with sine window */ - int i; - double xi[n]; - - for (i = 0; i < n; i++) { - xr[i] = sin((double)i/(double)n*M_PI)*xr[i]; xi[i] = 0.0; - } - FFT(n, xr, xi); - for (i = 0; i < n; i++) - xr[i] = sqrt(xr[i]*xr[i]+xi[i]*xi[i]); -} - -void sin_FFT(int n, double xr[]) -{ - /* DFT with sine window */ - int i; - double *xi; - - xi = dvector(1, 2*n); - - for (i = 1; i <= n; i++) { - xi[2*i-1] = sin((double)i/n*M_PI)*xr[i-1]; xi[2*i] = 0.0; - } - dfour1(xi, (unsigned long)n, 1); - for (i = 1; i <= n; i++) - xr[i-1] = sqrt(pow(xi[2*i-1], 2)+pow(xi[2*i], 2))*2.0/n; - - free_dvector(xi, 1, 2*n); -} - - -void sin_FFT(int n, double xr[], double xi[]) -{ - /* DFT with sine window */ - int i; - double *xri; - - xri = dvector(1, 2*n); - - for (i = 1; i <= n; i++) { - xri[2*i-1] = sin((double)i/n*M_PI)*xr[i-1]; - xri[2*i] = sin((double)i/n*M_PI)*xi[i-1]; - } - dfour1(xri, (unsigned long)n, 1); - for (i = 1; i <= n; i++) { - xr[i-1] = sqrt(pow(xri[2*i-1], 2)+pow(xri[2*i], 2))*2.0/n; - xi[i-1] = atan2(xri[2*i], xri[2*i-1]); - } - - free_dvector(xri, 1, 2*n); -} - - -void GetInd(int n, int k, int *ind1, int *ind3) -{ - if (k == 1) { - *ind1 = 2; *ind3 = 2; - } else if (k == n/2+1) { - *ind1 = n/2; *ind3 = n/2; - } else { - *ind1 = k - 1; *ind3 = k + 1; - } -} - - -void GetInd1(int n, int k, int *ind1, int *ind3) -{ - if (k == 1) { - *ind1 = 2; *ind3 = 2; - } else if (k == n) { - *ind1 = n - 1; *ind3 = n - 1; - } else { - *ind1 = k - 1; *ind3 = k + 1; - } -} - - -void GetPeak(int n, double *x, int *k) -{ - /* Locate peak in DFT spectrum */ - int ind1, ind2, ind3; - double peak; - - peak = 0.0; *k = 1; - for (ind2 = 1; ind2 <= n/2+1; ind2++) { - GetInd(n, ind2, &ind1, &ind3); - if (x[ind2-1] > peak && x[ind1-1] < x[ind2-1] && - x[ind3-1] < x[ind2-1]) - { - peak = x[ind2-1]; *k = ind2; - } - } -} - - -void GetPeak1(int n, double *x, int *k) -{ - /* Locate peak in DFT spectrum */ - int ind1, ind2, ind3; - double peak; - - peak = 0.0; *k = 1; - for (ind2 = 1; ind2 <= n; ind2++) { - GetInd1(n, ind2, &ind1, &ind3); - if (x[ind2-1] > peak && x[ind1-1] < x[ind2-1] && - x[ind3-1] < x[ind2-1]) - { - peak = x[ind2-1]; *k = ind2; - } - } -} - - -double Int2snu(int n, double *x, int k) -{ - /* Get frequency by nonlinear interpolation with two samples - for sine window. The interpolation is: - - 1 2 A(k) 1 - nu = - [ k - 1 + ------------- - - ] , k-1 <= N nu <= k - N A(k-1) + A(k) 2 - */ - int ind, ind1, ind3; - double ampl1, ampl2; - - GetInd(n, k, &ind1, &ind3); - if (x[ind3 - 1] > x[ind1 - 1]) { - ampl1 = x[k - 1]; - ampl2 = x[ind3 - 1]; - ind = k; - } else { - ampl1 = x[ind1 - 1]; - ampl2 = x[k - 1]; - /* Interpolate in right direction for 0 frequency */ - if (k != 1) - ind = ind1; - else - ind = 0; - } - /* Avoid division by zero */ - if (ampl1 + ampl2 != 0.0) - return ((ind - 1 + 2 * ampl2 / (ampl1 + ampl2) - 0.5) / n); - else - return 0.0; -} - - -/****************************************************************************/ -/* double Sinc(double omega) - - Purpose: - - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - zero test to be changed numericallywise - -****************************************************************************/ -double Sinc(double omega) -{ - /* Function to calculate: - - sin( omega ) - ------------ - omega - */ - if (omega != 0.0) - return (sin(omega) / omega); - else - return 1.0; -} - - -/****************************************************************************/ -/* double intsampl(long n, double *x, double nu, long k) - - Purpose: - - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - 17/07/03 use M_PI instead of pi - -****************************************************************************/ -double intsampl(int n, double *x, double nu, int k) -{ - /* Get amplitude by nonlinear interpolation for sine window. The - distribution is given by: - - 1 sin pi ( k + 1/2 ) sin pi ( k - 1/2 ) - F(k) = - ( -------------------- + -------------------- ) - 2 pi ( k + 1/2 ) pi ( k - 1/2 ) - */ - double corr; - - corr = (Sinc(M_PI * (k - 1 + 0.5 - nu * n)) + - Sinc(M_PI * (k - 1 - 0.5 - nu * n))) / 2; - return (x[k - 1] / corr); -} - - -/****************************************************************************/ -/* double linint(long n, long k, double nu, double *x) - - Purpose: - - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - 17/07/03 use M_PI instead of pi - -****************************************************************************/ -double linint(int n, int k, double nu, double *x) -{ - /* Get phase by linear interpolation for rectangular window - with -pi <= phi <= pi */ - int i; - double phi; - double xr[n], xi[n]; - - for (i = 0; i < n; i++) { - xr[i] = x[i]; xi[i] = 0.0; - } - FFT(n, xr, xi); - phi = GetAngle(xr[k-1], xi[k-1]) - (n*nu-k+1)*M_PI; - if (phi > M_PI) - phi -= 2.0*M_PI; - else if (phi < -M_PI) - phi += 2.0*M_PI; - return phi; -} - -/****************************************************************************/ -/* void FndRes(struct LOC_findres *LINK) - - Purpose: - - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - none - -****************************************************************************/ -void FndRes(struct LOC_findres *LINK) -{ - int i, j, FORLIM, FORLIM1; - double delta; - - FORLIM = LINK->n; - for (i = 0; i <= FORLIM; i++) { - FORLIM1 = LINK->n; - for (j = -LINK->n; j <= FORLIM1; j++) { - delta = fabs(i * LINK->nux + j * LINK->nuy); - delta -= (int)delta; - if (delta > 0.5) - delta = 1 - delta; - delta = fabs(delta - LINK->f); - delta -= (int)delta; - if (delta > 0.5) - delta = 1 - delta; - if (delta < LINK->eps) { - if (abs(i) + abs(j) < LINK->n && (i != 0 || j >= 0)) { - LINK->found = true; - *LINK->nx = i; - *LINK->ny = j; - } - } - } - } -} - - -/****************************************************************************/ -/* void FindRes(long n_, double nux_, double nuy_, double f_, - long *nx_, long *ny_) - - Purpose: - - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - none - -****************************************************************************/ -void FindRes(int n_, double nux_, double nuy_, double f_, int *nx_, int *ny_) -{ - /* Match f by a linear combination of nux and nuy */ - struct LOC_findres V; - - V.n = n_; - V.nux = nux_; - V.nuy = nuy_; - V.f = f_; - V.nx = nx_; - V.ny = ny_; - V.found = false; - V.eps = 0.5e-6; - do { - V.eps = 10 * V.eps; - FndRes(&V); - } while (!V.found); -} - - -void GetPeaks(int n, double *x, int nf, double *nu, double *A) -{ - int i, k, ind1, ind3; - - for (i = 0; i < nf; i++) { - GetPeak(n, x, &k); - nu[i] = Int2snu(n, x, k); A[i] = intsampl(n, x, nu[i], k); - /* Make peak flat to allow for new call */ - GetInd(n, k, &ind1, &ind3); - if (x[ind1-1] > x[ind3-1]) - x[k-1] = x[ind1-1]; - else - x[k-1] = x[ind3-1]; - } -} - - -void GetPeaks1(int n, double *x, int nf, double *nu, double *A) -{ - int i, k, ind1, ind3; - - for (i = 0; i < nf; i++) { - GetPeak1(n, x, &k); - nu[i] = Int2snu(n, x, k); A[i] = intsampl(n, x, nu[i], k); - /* Make peak flat to allow for new call */ - GetInd1(n, k, &ind1, &ind3); - if (x[ind1-1] > x[ind3-1]) - x[k-1] = x[ind1-1]; - else - x[k-1] = x[ind3-1]; - } -} - -/*******************************/ -/* Routines for magnetic error */ -/*******************************/ - -void SetTol(int Fnum, double dxrms, double dyrms, double drrms) -{ - int i; - long k; - - for (i = 1; i <= GetnKid(Fnum); i++) { - k = Elem_GetPos(Fnum, i); - Cell[k].Elem.M->PdSrms[X_] = dxrms; - Cell[k].Elem.M->PdSrnd[X_] = normranf(); - Cell[k].Elem.M->PdSrms[Y_] = dyrms; - Cell[k].Elem.M->PdSrnd[Y_] = normranf(); - Cell[k].Elem.M->PdTrms = drrms; - Cell[k].Elem.M->PdTrnd = normranf(); - Mpole_SetdS(Fnum, i); Mpole_SetdT(Fnum, i); - } -} - - -void Scale_Tol(int Fnum, double dxrms, double dyrms, double drrms) -{ - int Knum; - long int loc; - - for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) { - loc = Elem_GetPos(Fnum, Knum); - Cell[loc].Elem.M->PdSrms[X_] = dxrms; Cell[loc].Elem.M->PdSrms[Y_] = dyrms; - Cell[loc].Elem.M->PdTrms = drrms; - Mpole_SetdS(Fnum, Knum); Mpole_SetdT(Fnum, Knum); - } -} - - -/****************************************************************************/ -/* void SetaTol(int Fnum, int Knum, double dx, double dy, double dr) - - Purpose: - Set a known random multipole displacement error - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - none - -****************************************************************************/ -void SetaTol(int Fnum, int Knum, double dx, double dy, double dr) -{ - long int loc; - - loc = Elem_GetPos(Fnum, Knum); - Cell[loc].Elem.M->PdSrms[0] = dx; Cell[loc].Elem.M->PdSrnd[0] = 1e0; - Cell[loc].Elem.M->PdSrms[1] = dy; Cell[loc].Elem.M->PdSrnd[1] = 1e0; - Cell[loc].Elem.M->PdTrms = dr; Cell[loc].Elem.M->PdTrnd = 1e0; - Mpole_SetdS(Fnum, Knum); Mpole_SetdT(Fnum, Knum); -} - - -void ini_aper(const double Dxmin, const double Dxmax, - const double Dymin, const double Dymax) -{ - int k; - - for (k = 0; k <= globval.Cell_nLoc; k++) { - Cell[k].maxampl[X_][0] = Dxmin; Cell[k].maxampl[X_][1] = Dxmax; - Cell[k].maxampl[Y_][0] = Dymin; Cell[k].maxampl[Y_][1] = Dymax; - } -} - -void set_aper(const int Fnum, const double Dxmin, const double Dxmax, - const double Dymin, const double Dymax) -{ - int i; - long int loc; - - for (i = 1; i <= GetnKid(Fnum); i++) { - loc = Elem_GetPos(Fnum, i); - Cell[loc].maxampl[X_][0] = Dxmin; Cell[loc].maxampl[X_][1] = Dxmax; - Cell[loc].maxampl[Y_][0] = Dymin; Cell[loc].maxampl[Y_][1] = Dymax; - } -} - - -void LoadApertures(const char *ChamberFileName) -{ - char line[128], FamName[32]; - long Fnum; - double Xmin, Xmax, Ymin, Ymax; - FILE *ChamberFile; - - ChamberFile = file_read(ChamberFileName); - - do - fgets(line, 128, ChamberFile); - while (strstr(line, "#") != NULL); - - do { - sscanf(line,"%s %lf %lf %lf %lf", FamName,&Xmin, &Xmax, &Ymin,&Ymax); - Fnum = ElemIndex(FamName); - if (Fnum > 0) set_aper(Fnum, Xmin, Xmax, Ymin, Ymax); - } while (fgets(line, 128, ChamberFile ) != NULL); - - fclose(ChamberFile); -} - - -// Load tolerances from the file -void LoadTolerances(const char *TolFileName) -{ - char line[128], FamName[32]; - int Fnum; - double dx, dy, dr; - FILE *tolfile; - - tolfile = file_read(TolFileName); - - do - fgets(line, 128, tolfile); - while (strstr(line, "#") != NULL); - - do { - if (strstr(line, "#") == NULL) { - sscanf(line,"%s %lf %lf %lf", FamName, &dx, &dy, &dr); - Fnum = ElemIndex(FamName); - if (Fnum > 0) { - SetTol(Fnum, dx, dy, dr); - } else { - printf("LoadTolerances: undefined element %s\n", FamName); - exit_(1); - } - } - } while (fgets(line, 128, tolfile) != NULL); - - fclose(tolfile); -} - - -// Load tolerances from the file -void ScaleTolerances(const char *TolFileName, const double scl) -{ - char line[128], FamName[32]; - int Fnum; - double dx, dy, dr; - FILE *tolfile; - - tolfile = file_read(TolFileName); - - do - fgets(line, 128, tolfile); - while (strstr(line, "#") != NULL); - - do { - if (strstr(line, "#") == NULL) { - sscanf(line,"%s %lf %lf %lf", FamName, &dx, &dy, &dr); - Fnum = ElemIndex(FamName); - if (Fnum > 0) { - Scale_Tol(Fnum, scl*dx, scl*dy, scl*dr); - } else { - printf("ScaleTolerances: undefined element %s\n", FamName); - exit_(1); - } - } - } while (fgets(line, 128, tolfile) != NULL); - fclose(tolfile); -} - - -void SetKpar(int Fnum, int Knum, int Order, double k) -{ - - Cell[Elem_GetPos(Fnum, Knum)].Elem.M->PBpar[Order+HOMmax] = k; - Mpole_SetPB(Fnum, Knum, Order); -} - - -void SetL(int Fnum, int Knum, double L) -{ - - Cell[Elem_GetPos(Fnum, Knum)].Elem.PL = L; -} - - -void SetL(int Fnum, double L) -{ - int i; - - for (i = 1; i <= GetnKid(Fnum); i++) - Cell[Elem_GetPos(Fnum, i)].Elem.PL = L; -} - - -void SetdKpar(int Fnum, int Knum, int Order, double dk) -{ - - Cell[Elem_GetPos(Fnum, Knum)].Elem.M->PBpar[Order+HOMmax] += dk; - Mpole_SetPB(Fnum, Knum, Order); -} - - -void SetKLpar(int Fnum, int Knum, int Order, double kL) -{ - long int loc; - - loc = Elem_GetPos(Fnum, Knum); - if (Cell[loc].Elem.PL != 0e0) - Cell[loc].Elem.M->PBpar[Order+HOMmax] = kL/Cell[loc].Elem.PL; - else - Cell[loc].Elem.M->PBpar[Order+HOMmax] = kL; - Mpole_SetPB(Fnum, Knum, Order); -} - - -void SetdKLpar(int Fnum, int Knum, int Order, double dkL) -{ - long int loc; - - loc = Elem_GetPos(Fnum, Knum); - if (Cell[loc].Elem.PL != 0e0) - Cell[loc].Elem.M->PBpar[Order + HOMmax] += dkL/Cell[loc].Elem.PL; - else - Cell[loc].Elem.M->PBpar[Order + HOMmax] += dkL; - Mpole_SetPB(Fnum, Knum, Order); -} - - -void SetdKrpar(int Fnum, int Knum, int Order, double dkrel) -{ - long int loc; - - loc = Elem_GetPos(Fnum, Knum); - if (Order == Dip && Cell[loc].Elem.M->Pthick == thick) - Cell[loc].Elem.M->PBpar[Dip+HOMmax] += dkrel*Cell[loc].Elem.M->Pirho; - else - Cell[loc].Elem.M->PBpar[Order+HOMmax] - += dkrel*Cell[loc].Elem.M->PBpar[Order+HOMmax]; - Mpole_SetPB(Fnum, Knum, Order); -} - - -void Setbn(int Fnum, int order, double bn) -{ - int i; - - for (i = 1; i <= GetnKid(Fnum); i++) - SetKpar(Fnum, i, order, bn); -} - - -void SetbnL(int Fnum, int order, double bnL) -{ - int i; - - for (i = 1; i <= GetnKid(Fnum); i++) - SetKLpar(Fnum, i, order, bnL); -} - - -void Setdbn(int Fnum, int order, double dbn) -{ - int i; - - for (i = 1; i <= GetnKid(Fnum); i++) - SetdKpar(Fnum, i, order, dbn); -} - - -void SetdbnL(int Fnum, int order, double dbnL) -{ - int i; - - for (i = 1; i <= GetnKid(Fnum); i++) { - SetdKLpar(Fnum, i, order, dbnL); - } -} - - -void Setbnr(int Fnum, long order, double bnr) -{ - int i; - - for (i = 1; i <= GetnKid(Fnum); i++) - SetdKrpar(Fnum, i, order, bnr); -} - - -void SetbnL_sys(int Fnum, int Order, double bnL_sys) -{ - int Knum; - long int loc; - - for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) { - loc = Elem_GetPos(Fnum, Knum); - if (Cell[loc].Elem.PL != 0.0) - Cell[loc].Elem.M->PBsys[Order+HOMmax] = bnL_sys/Cell[loc].Elem.PL; - else - Cell[loc].Elem.M->PBsys[Order+HOMmax] = bnL_sys; - Mpole_SetPB(Fnum, Knum, Order); - } -} - - -void set_dbn_rel(const int type, const int n, const double dbn_rel) -{ - long int j; - double dbn; - - printf("\n"); - printf("Setting Db_%d/b_%d = %6.1e for:\n", n, type, dbn_rel); - printf("\n"); - for (j = 0; j <= globval.Cell_nLoc; j++) - if ((Cell[j].Elem.Pkind == Mpole) && (Cell[j].Elem.M->n_design == type)) { - printf("%s\n", Cell[j].Elem.PName); - dbn = dbn_rel*Cell[j].Elem.M->PBpar[type+HOMmax]; - Cell[j].Elem.M->PBrms[n+HOMmax] = dbn; - Cell[j].Elem.M->PBrnd[n+HOMmax] = normranf(); - Mpole_SetPB(Cell[j].Fnum, Cell[j].Knum, n); - } -} - - -double GetKpar(int Fnum, int Knum, int Order) -{ - return (Cell[Elem_GetPos(Fnum, Knum)].Elem.M->PBpar[Order+HOMmax]); -} - - -double GetL(int Fnum, int Knum) -{ - return (Cell[Elem_GetPos(Fnum, Knum)].Elem.PL); -} - - -double GetKLpar(int Fnum, int Knum, int Order) -{ - long int loc; - - loc = Elem_GetPos(Fnum, Knum); - if (Cell[loc].Elem.PL != 0e0) - return (Cell[loc].Elem.M->PBpar[Order+HOMmax]*Cell[loc].Elem.PL); - else - return (Cell[loc].Elem.M->PBpar[Order+HOMmax]); -} - - -void SetdKLrms(int Fnum, int Order, double dkLrms) -{ - long int Knum, loc; - - for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) { - loc = Elem_GetPos(Fnum, Knum); - if (Cell[loc].Elem.PL != 0e0) - Cell[loc].Elem.M->PBrms[Order+HOMmax] = dkLrms/Cell[loc].Elem.PL; - else - Cell[loc].Elem.M->PBrms[Order+HOMmax] = dkLrms; - Cell[loc].Elem.M->PBrnd[Order+HOMmax] = normranf(); - Mpole_SetPB(Fnum, Knum, Order); - } -} - - -void Setdkrrms(int Fnum, int Order, double dkrrms) -{ - long int Knum, loc; - - for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) { - loc = Elem_GetPos(Fnum, Knum); - if (Order == Dip && Cell[loc].Elem.M->Pthick == thick) - Cell[loc].Elem.M->PBrms[Dip+HOMmax] = dkrrms*Cell[loc].Elem.M->Pirho; - else - Cell[loc].Elem.M->PBrms[Order+HOMmax] - = dkrrms*Cell[loc].Elem.M->PBpar[Order+HOMmax]; - Cell[loc].Elem.M->PBrnd[Order+HOMmax] = normranf(); - Mpole_SetPB(Fnum, Knum, Order); - } -} - - -void SetKL(int Fnum, int Order) -{ - long int Knum; - - for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) - Mpole_SetPB(Fnum, Knum, Order); -} - - -void set_dx(const int type, const double sigma_x, const double sigma_y) -{ - long int j; - - printf("\n"); - printf("Setting sigma_x,y = (%6.1e, %6.1e) for b_%d:\n", - sigma_x, sigma_y, type); - printf("\n"); - for (j = 0; j <= globval.Cell_nLoc; j++) - if ((Cell[j].Elem.Pkind == Mpole) && (Cell[j].Elem.M->n_design == type)) { - printf("%s\n", Cell[j].Elem.PName); - Cell[j].Elem.M->PdSrms[X_] = sigma_x; - Cell[j].Elem.M->PdSrms[Y_] = sigma_y; - Cell[j].Elem.M->PdSrnd[X_] = normranf(); - Cell[j].Elem.M->PdSrnd[Y_] = normranf(); - Mpole_SetdS(Cell[j].Fnum, Cell[j].Knum); - } -} - - -void SetBpmdS(int Fnum, double dxrms, double dyrms) -{ - long int Knum, loc; - - for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) { - loc = Elem_GetPos(Fnum, Knum); - Cell[loc].dS[X_] = normranf()*dxrms; Cell[loc].dS[Y_] = normranf()*dyrms; - } -} - - -/****************************************/ -/* Routines for closed orbit correction */ -/****************************************/ - - -/****************************************************************************/ -void codstat(double *mean, double *sigma, double *xmax, long lastpos, bool all) -{ - long i, j, n; - Vector2 sum, sum2; - double TEMP; - - n = 0; - for (j = 0; j <= 1; j++) { - sum[j] = 0.0; sum2[j] = 0.0; xmax[j] = 0.0; - } - for (i = 0; i <= lastpos; i++) { - if (all || Cell[i].Fnum == globval.bpm) { - n++; - for (j = 1; j <= 2; j++) { - sum[j - 1] += Cell[i].BeamPos[j * 2 - 2]; - TEMP = Cell[i].BeamPos[j * 2 - 2]; - sum2[j - 1] += TEMP * TEMP; - xmax[j - 1] = max(xmax[j - 1], fabs(Cell[i].BeamPos[j * 2 - 2])); - } - } - } - for (j = 0; j <= 1; j++) { - if (n != 0) - mean[j] = sum[j] / n; - else - mean[j] = 0.0; - if (n != 0 && n != 1) { - TEMP = sum[j]; - sigma[j] = (n * sum2[j] - TEMP * TEMP) / (n * (n - 1.0)); - } else - sigma[j] = 0.0; - if (sigma[j] >= 0.0) - sigma[j] = sqrt(sigma[j]); - else - sigma[j] = 0.0; - } -} - -/****************************************************************************/ -/* void CodStatBpm(double *mean, double *sigma, double *xmax, long lastpos, - long bpmdis[mnp]) - - Purpose: - Get statistics for closed orbit - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - none - -****************************************************************************/ -void CodStatBpm(double *mean, double *sigma, double *xmax, long lastpos, - long bpmdis[mnp]) -{ - long i, j, m, n; - Vector2 sum, sum2; - double TEMP; - - m= n= 0; - for (j = 0; j <= 1; j++) { - sum[j] = 0.0; sum2[j] = 0.0; xmax[j] = 0.0; - } - - for (i = 0; i <= lastpos; i++) { - if (Cell[i].Fnum == globval.bpm) { - if (! bpmdis[m]) { - for (j = 1; j <= 2; j++) { - sum[j - 1] += Cell[i].BeamPos[j * 2 - 2]; - TEMP = Cell[i].BeamPos[j * 2 - 2]; - sum2[j - 1] += TEMP * TEMP; - xmax[j - 1] = max(xmax[j - 1], fabs(Cell[i].BeamPos[j * 2 - 2])); - } - n++; - } - m++; - } - } - for (j = 0; j <= 1; j++) { - if (n != 0) - mean[j] = sum[j] / n; - else - mean[j] = 0.0; - if (n != 0 && n != 1) { - TEMP = sum[j]; - sigma[j] = (n * sum2[j] - TEMP * TEMP) / (n * (n - 1.0)); - } else - sigma[j] = 0.0; - if (sigma[j] >= 0.0) - sigma[j] = sqrt(sigma[j]); - else - sigma[j] = 0.0; - } -} - - - -/****************************************************************************/ -/* double digitize(double x, double maxkick, double maxsamp) - - Purpose: - Map x onto the integer interval (-maxsamp ... maxsamp) where maxsamp - corresponds maxkick. - - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - none - -****************************************************************************/ -double digitize(double x, double maxkick, double maxsamp) -{ - if (maxkick>0.) - if (maxsamp>1.) - return Sgn(x)*maxkick/maxsamp *min(floor(fabs(x)/maxkick*maxsamp),maxsamp-1.); - else { - return Sgn(x)*min(fabs(x),maxkick); - } - else - return x; -} - - -/****************************************************************************/ -/* double digitize2(long plane, long inum, double x, double maxkick, double maxsamp) - - Purpose: - - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - none - -****************************************************************************/ -svdarray xmemo[2]; - -double digitize2(long plane, long inum, double x, double maxkick, double maxsamp) -{ - double xint; - - if (maxkick>0.) - if (maxsamp>1.) - { - xint=min(floor(fabs(x)/maxkick*maxsamp),maxsamp-1.); - - if(fabs(xint-xmemo[inum][plane]) >=1.) - { - xmemo[inum][plane]=xint; - } - else - { - xmemo[inum][plane]+=0.1; - xint=min(xmemo[inum][plane],maxsamp-1.); - } - return Sgn(x)*maxkick/maxsamp*xint; - } - else - { - return Sgn(x)*min(fabs(x),maxkick); - } - else - return x; -} - - -// MATH ROUTINE a mettre dans mathlib.c - -/****************************************************************************/ -/* void GetMean(n, x) - - Purpose: - Get out the mean value of vector x - - Input: - n vector size - x vector to get out the mean value - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - to be moved in mathlib - -****************************************************************************/ -void GetMean(long n, double *x) -{ - long i; - double mean = 0e0; - - if ( n < 1 ) - { - fprintf(stdout,"GetMean: error wrong vector size n=%ld\n",n); - exit_(1); - } - for (i = 0; i < n; i++) - mean += x[i]; - mean /= n; - for (i = 0; i < n; i++) - x[i] = x[i] - mean; -} - -/****************************************************************************/ -/* double Fract(double x) - - Purpose: - Gets fractional part of x - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - none - -****************************************************************************/ -double Fract(double x) -{ - return (x - (long int)x); -} - -/****************************************************************************/ -/* double Sgn (double x) - - Purpose: - Gets sign of x - - Input: - none - - Output: - 0 if zero - 1 if positive - -1 if negative - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - none - -****************************************************************************/ -double Sgn (double x) -{ - return (x == 0.0 ? 0.0 : (x > 0.0 ? 1.0 : -1.0)); -} - - -void ChamberOff(void) -{ - int i; - - for (i = 0; i <= globval.Cell_nLoc; i++) { - Cell[i].maxampl[X_][0] = -max_ampl; Cell[i].maxampl[X_][1] = max_ampl; - Cell[i].maxampl[Y_][0] = -max_ampl; Cell[i].maxampl[Y_][1] = max_ampl; - } - status.chambre = false; -} - - -void PrintCh(void) -{ - long i = 0; - struct tm *newtime; - FILE *f; - - const char *fic = "chambre.out"; - - newtime = GetTime(); - - f = file_write(fic); - fprintf(f, "# TRACY II v.2.6 -- %s -- %s \n", fic, asctime2(newtime)); - fprintf(f, "# name s -xch +xch zch\n"); - fprintf(f, "# [mm] [mm] [mm]\n"); - fprintf(f, "#\n"); - - for (i = 0; i <= globval.Cell_nLoc; i++) - fprintf(f, "%4ld %15s %6.2f %7.3f %7.3f %7.3f\n", - i, Cell[i].Elem.PName, Cell[i].S, - Cell[i].maxampl[X_][0]*1E3, Cell[i].maxampl[X_][1]*1E3, - Cell[i].maxampl[Y_][1]*1E3); - - fclose(f); -} - - -/** function from soleilcommon.c **/ - -void Read_Lattice(char *fic) -{ - bool status; - char fic_maille[S_SIZE+4] = "", fic_erreur[S_SIZE+4] = ""; - int i; - double dP = 0.0; - Vector2 beta, alpha, eta, etap; - Vector codvect; - - const double RFacceptance = 0.060001; // soleil energy acceptance - - strcpy(fic_maille, fic); strcpy(fic_erreur, fic); - - /* generation automatique du nom du fichier maille et erreur */ - strcat(fic_maille, ".lat"); strcat(fic_erreur, ".lax"); - - /* Initialisation de Tracy */ - - t2init(); - - /* open the lattice Input file */ - - if ((fi = fopen(fic_maille, "r")) == NULL) { - fprintf(stdout, "ReadLattice: Error while opening file %s \n", fic_maille); - fprintf(stdout, "The lattice file has to end by .lat \n"); - exit_(1); - } - - /* opens the lattice Output file */ - if ((fo = fopen(fic_erreur, "w")) == NULL) { - fprintf(stdout, "ReadLattice: Error while opening file %s \n", fic_erreur); - exit_(1); - } - - /* Reads lattice and set principle parameters - * Energy CODeps and energy offset - * print statistics - */ - status = Lattice_Read(&fi, &fo); - - if (status == false) { - cout << "Lattice_Read function has returned false" << endl; - cout << "See file " << fic_erreur << endl; - exit_(1); - } - cout << "Lattice file: " << fic_maille << endl; - - /* initializes cell structure: construction of the RING */ - /* Creator of all the matrices for each element */ - Cell_Init(); - - if (globval.RingType == 1) { // for a ring - /* define x/y physical aperture */ - ChamberOff(); - - /* Defines global variables for Tracy code */ - globval.H_exact = false; // Small Ring Hamiltonian - globval.quad_fringe = false; // quadrupole fringe fields on/off - globval.EPU = false; // Elliptically Polarizing Undulator - globval.Cavity_on = false; /* Cavity on/off */ - globval.radiation = false; /* radiation on/off */ - globval.IBS = false; /* diffusion on/off */ - globval.emittance = false; /* emittance on/off */ - globval.pathlength = false; /* Path lengthening computation */ - globval.CODimax = 40; /* maximum number of iterations for COD - algo */ - globval.delta_RF = RFacceptance;/* energy acceptance for SOLEIL */ - - Cell_SetdP(dP); /* added for correcting BUG if non convergence: - compute on momentum linear matrices */ - } else { - // for transfer lines - /* Initial settings : */ - beta[X_] = 8.1; alpha[X_] = 0.0; beta[Y_] = 8.1; alpha[Y_] = 0.0; - eta[X_] = 0.0; etap[X_] = 0.0; eta[Y_] = 0.0; etap[Y_] = 0.0; - - for (i = 0; i < ss_dim; i++) { - codvect[i] = 0.0; globval.CODvect[i] = codvect[i]; - } - dP = codvect[delta_]; - - /* Defines global variables for Tracy code */ - globval.Cavity_on = false; /* Cavity on/off */ - globval.radiation = false; /* radiation on/off */ - globval.emittance = false; /* emittance on/off */ - globval.pathlength = false; /* Path lengthening computation */ - globval.CODimax = 10; /* maximum number of iterations for COD - algo */ - globval.delta_RF = RFacceptance; /* 6% + epsilon energy acceptance - for SOLEIL */ - globval.dPparticle = dP; - - ChamberOff(); - - TransTwiss(alpha, beta, eta, etap, codvect); - } -} - - -/****************************************************************************/ -/* void GetChromTrac(long Nb, long Nbtour, double emax, - double *xix, double *xiz) - - Purpose: - Computes chromaticities by tracking - - Input: - Nb point number - Nbtour turn number - emax energy step - - Output: - xix horizontal chromaticity - xiz vertical chromaticity - - Return: - none - - Global variables: - trace - - Specific functions: - Trac_Simple, Get_NAFF - - Comments: - 27/04/03 chromaticities are now output arguments - -****************************************************************************/ -#define nterm 2 -void GetChromTrac(long Nb, long Nbtour, double emax, double *xix, double *xiz) -{ - bool status = true; - int nb_freq[2] = { 0, 0 }; /* frequency number to look for */ - int i = 0; - double Tab[6][NTURN], fx[nterm], fz[nterm], nux1, nux2, nuz1, nuz2; - - double x = 1e-6, xp = 0.0, z = 1e-6, zp = 0.0; - double x0 = 1e-6, xp0 = 0.0, z0 = 1e-6, zp0 = 0.0; - - /* initializations */ - for (i = 0; i < nterm; i++) { - fx[i] = 0.0; fz[i] = 0.0; - } - /* end init */ - - /* Tracking for delta = emax and computing tunes */ - x = x0; xp = xp0; z = z0; zp = zp0; - - Trac_Simple(x, xp, z, zp, emax, 0.0, Nbtour, Tab, &status); - Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq); - - nux1 = (fabs (fx[0]) > 1e-8 ? fx[0] : fx[1]); nuz1 = fz[0]; - - if (trace) - fprintf(stdout, - "\n Entering routine for chroma using tracking\n"); - if (trace) - fprintf(stdout, "emax= % 10.6e nux1=% 10.6e nuz1= % 10.6e\n", - emax, nux1, nuz1); - - /* Tracking for delta = -emax and computing tunes */ - x = x0; xp = xp0; z = z0; zp = zp0; - - Trac_Simple(x, xp, z, zp, -emax, 0.0, Nbtour, Tab, &status); - Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq); - - if (trace) - fprintf(stdout, "nturn=%6ld x=% 10.5g xp=% 10.5g z=% 10.5g zp=% 10.5g" - " delta=% 10.5g ctau=% 10.5g \n", - Nbtour, - Tab[0][Nbtour-1], Tab[1][Nbtour-1], - Tab[2][Nbtour-1], Tab[3][Nbtour-1], - Tab[4][Nbtour-1], Tab[5][Nbtour-1]); - - nux2 = (fabs(fx[0]) > 1e-8 ? fx[0] : fx[1]); nuz2 = fz[0]; - - if (trace) - fprintf(stdout, "emax= % 10.6e nux2= % 10.6e nuz2= % 10.6e\n", - -emax, nux2, nuz2); - - /* Computing chromaticities */ - *xix = (nux2-nux1)*0.5/emax; *xiz = (nuz2-nuz1)*0.5/emax; - - if (trace) - fprintf (stdout, " Exiting routine for chroma using tracking\n\n"); -} -#undef nterm - -/****************************************************************************/ -/* void GetTuneTrac(long Nbtour, double emax, double *nux, double *nuz) - - Purpose: - Computes chromaticities by tracking - - Input: - Nb point number - Nbtour turn number - emax energy step - - Output: - none - - Return: - none - - Global variables: - trace - - Specific functions: - Trac_Simple, Get_NAFF - - Comments: - none - -****************************************************************************/ -#define nterm 2 -void GetTuneTrac(long Nbtour, double emax, double *nux, double *nuz) -{ - double Tab[6][NTURN], fx[nterm], fz[nterm]; - int nb_freq[2]; - bool status; - - double x = 1e-6, xp = 0.0, z = 1e-6, zp = 0.0; - - Trac_Simple(x, xp, z, zp, emax, 0.0, Nbtour, Tab, &status); - Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq); - - *nux = (fabs (fx[0]) > 1e-8 ? fx[0] : fx[1]); - *nuz = fz[0]; -} -#undef nterm - -/****************************************************************************/ -/* void TransTwiss(double *alpha, double *beta, double *eta, double *etap, double *codvect) - - Purpose: high level application - Calculate Twiss functions for a transport line - - Input: - alpha alpha fonctions at the line entrance - beta beta fonctions at the line entrance - eta disperion fonctions at the line entrance - etap dipersion derivatives fonctions at the line entrance - codvect closed orbit fonctions at the line entrance - - Output: - none - - Return: - none - - Global variables: - - - Specific functions: - TransTrace - - Comments: - redundant with ttwiss - -****************************************************************************/ -void TransTwiss(Vector2 &alpha, Vector2 &beta, Vector2 &eta, Vector2 &etap, - Vector &codvect) -{ - TransTrace(0, globval.Cell_nLoc, alpha, beta, eta, etap, codvect); -} - - -/****************************************************************************/ -/* void ttwiss(double *alpha, double *beta, double *eta, double *etap, double dP) - - Purpose: - Calculate Twiss functions for transport line - - Input: - none - - Output: - none - - Return: - none - - Global variables: - none - - Specific functions: - none - - Comments: - redundant with TransTwiss - -****************************************************************************/ -void ttwiss(const Vector2 &alpha, const Vector2 &beta, - const Vector2 &eta, const Vector2 &etap, const double dP) -{ - TraceABN(0, globval.Cell_nLoc, alpha, beta, eta, etap, dP); -} - -/****************************************************************************/ -/* void findcodS(double dP) - - Purpose: - Search for the closed orbit using a numerical method - Algo: Newton_Raphson method - Quadratic convergence - May need a guess starting point - Simple precision algorithm - - Input: - dP energy offset - - Output: - none - - Return: - none - - Global variables: - none - - specific functions: - Newton_Raphson - - Comments: - Method introduced because of bad convergence of da for ID using RADIA maps - -****************************************************************************/ -void findcodS(double dP) -{ - double *vcod; - Vector x0; - const int ntrial = 40; // maximum number of trials for closed orbit - const double tolx = 1e-8; // numerical precision - int k; - int dim; // 4D or 6D tracking - long lastpos; - - vcod = dvector(1, 6); - - // starting point - for (k = 1; k <= 6; k++) - vcod[k] = 0.0; - - vcod[5] = dP; // energy offset - - if (globval.Cavity_on){ - dim = 6; /* 6D tracking*/ - fprintf(stdout,"Error looking for cod in 6D\n"); - exit_(1); - } - else{ - dim = 4; /* 4D tracking */ - vcod[1] = Cell[0].Eta[0]*dP; vcod[2] = Cell[0].Etap[0]*dP; - vcod[3] = Cell[0].Eta[1]*dP; vcod[4] = Cell[0].Etap[1]*dP; - } - - Newton_RaphsonS(ntrial, vcod, dim, tolx); - - if (status.codflag == false) - fprintf(stdout, "Error No COD found\n"); - if (trace) { - for (k = 1; k <= 6; k++) - x0[k-1] = vcod[k]; - fprintf(stdout, "Before cod % .5e % .5e % .5e % .5e % .5e % .5e \n", - x0[0], x0[1], x0[2], x0[3], x0[4], x0[5]); - Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); - fprintf(stdout, "After cod % .5e % .5e % .5e % .5e % .5e % .5e \n", - x0[0], x0[1], x0[2], x0[3], x0[4], x0[5]); - Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); - } - free_dvector(vcod,1,6); -} - -/****************************************************************************/ -/* void findcod(double dP) - - Purpose: - Search for the closed orbit using a numerical method - Algo: Newton_Raphson method - Quadratic convergence - May need a guess starting point - Simple precision algorithm - 4D - Starting point: linear closed orbit - - 6D - Starting point: zero - if radiation on : x[5] is the synchroneous phase (equilibrium RF phase) - off: x[5] is zero - - Input: - dP energy offset - - Output: - none - - Return: - none - - Global variables: - none - - specific functions: - Newton_Raphson - - Comments: - Method introduced because of bad convergence of da for ID - using RADIA maps - -****************************************************************************/ -void findcod(double dP) -{ - Vector vcod; - const int ntrial = 40; // maximum number of trials for closed orbit - const double tolx = 1e-10; // numerical precision - int k, dim = 0; - long lastpos; - - // initializations - for (k = 0; k <= 5; k++) - vcod[k] = 0.0; - - if (globval.Cavity_on){ - fprintf(stdout,"warning looking for cod in 6D\n"); - dim = 6; - } else{ // starting point linear closed orbit - dim = 4; - vcod[0] = Cell[0].Eta[0]*dP; vcod[1] = Cell[0].Etap[0]*dP; - vcod[2] = Cell[0].Eta[1]*dP; vcod[3] = Cell[0].Etap[1]*dP; - vcod[4] = dP; // energy offset - } - - Newton_Raphson(dim, vcod, ntrial, tolx); - - if (status.codflag == false) - fprintf(stdout, "Error No COD found\n"); - - CopyVec(6, vcod, globval.CODvect); // save closed orbit at the ring entrance - - if (trace) - { - fprintf(stdout, - "Before cod2 % .5e % .5e % .5e % .5e % .5e % .5e \n", - vcod[0], vcod[1], vcod[2], vcod[3], vcod[4], vcod[5]); - Cell_Pass(0, globval.Cell_nLoc, vcod, lastpos); - fprintf(stdout, - "After cod2 % .5e % .5e % .5e % .5e % .5e % .5e \n", - vcod[0], vcod[1], vcod[2], vcod[3], vcod[4], vcod[5]); - } -} -/****************************************************************************/ -/* void computeFandJS(double *x, int n, double **fjac, double *fvect) - - Purpose: - Simple precision algo - Tracks x over one turn. And computes the Jacobian matrix of the - transformation by numerical differentiation. - using forward difference formula : faster but less accurate - using symmetric difference formula - - Input: - x vector for evaluation - n dimension 4 or 6 - - Output: - fvect transport of x over one turn - fjac Associated jacobian matrix - - Return: - none - - Global variables: - none - - specific functions: - none - - Comments: - none - -****************************************************************************/ - -void computeFandJS(double *x, int n, double **fjac, double *fvect) -{ - int i, k; - long lastpos = 0L; - Vector x0, fx, fx1, fx2; - - const double deps = 1e-8; //stepsize for numerical differentiation - - for (i = 1; i <= 6; i++) - x0[i - 1] = x[i]; - - Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); - - for (i = 1; i <= n; i++) - { - fvect[i] = x0[i - 1]; - fx[i - 1] = x0[i - 1]; - } - - // compute Jacobian matrix by numerical differentiation - for (k = 0; k < n; k++) - { - for (i = 1; i <= 6; i++) - x0[i - 1] = x[i]; - x0[k] += deps; // differential step in coordinate k - - Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); // tracking along the ring - for (i = 1; i <= 6; i++) - fx1[i - 1] = x0[i - 1]; - - for (i = 1; i <= 6; i++) - x0[i - 1] = x[i]; - x0[5] = 0.0; - x0[k] -= deps; // differential step in coordinate k - - Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); // tracking along the ring - for (i = 1; i <= 6; i++) - fx2[i - 1] = x0[i - 1]; - - for (i = 1; i <= n; i++) // symmetric difference formula - fjac[i][k + 1] = 0.5 * (fx1[i - 1] - fx2[i - 1]) / deps; - //~ for (i = 1; i <= n; i++) // forward difference formula - //~ fjac[i][k + 1] = (float) ((x0[i - 1] - fx[i - 1]) / deps); - } -} - -/****************************************************************************/ -/* void computeFand(int n, float *x, float **fjac, float *fvect) - - Purpose: - Tracks x over one turn. And computes the Jacobian matrix of the - transformation by numerical differentiation. - using symmetric difference formula - double precision algorithm - - Input: - x vector for evaluation - - Output: - fvect transport of x over one turn - fjac Associated jacobian matrix - - Return: - none - - Global variables: - none - - specific functions: - none - - Comments: - none - -****************************************************************************/ -void computeFandJ(int n, Vector &x, Matrix &fjac, Vector &fvect) -{ - int i, k; - long lastpos = 0; - Vector x0, fx1, fx2; - - const double deps = 1e-8; //stepsize for numerical differentiation - - CopyVec(6, x, x0); - - Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); - CopyVec(n, x0, fvect); - - // compute Jacobian matrix by numerical differentiation - for (k = 0; k < n; k++) { - CopyVec(6L, x, x0); - x0[k] += deps; // differential step in coordinate k - - Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); // tracking along the ring - CopyVec(6L, x0, fx1); - - CopyVec(6L, x, x0); - x0[k] -= deps; // differential step in coordinate k - - Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); // tracking along the ring - CopyVec(6L, x0, fx2); - - for (i = 0; i < n; i++) // symmetric difference formula - fjac[i][k] = 0.5 * (fx1[i] - fx2[i]) / deps; - } -} - -/****************************************************************************/ -/* void Newton_RaphsonS(int ntrial,double x[],int n,double tolx, double tolf) - - Purpose: - Newton_Rapson algorithm from Numerical Recipes - single precision algorithm - Robustess: quadratic convergence - Hint: for n-dimensional problem, the algo can be stuck on local minimum - In this case, it should be enough to provide a resonable starting - point. - - Method: - look for closed orbit solution of f(x) = x - This problems is equivalent to finding the zero of g(x)= f(x) - x - g(x+h) ~= f(x) - x + (Jacobian(f) -Id) h + O(h*h) - Then at first order we solve h: - h = - inverse(Jacobian(f) -Id) * (f(x)-x) - the new guess is then xnew = x + h - By iteration, this converges quadratically. - - The algo is stopped whenever |x -xnew| < tolx - - f(x) is computes by tracking over one turn - Jacobian(f) is computed numerically by numerical differentiation - These two operations are provided by the function computeFandJ - - Input: - ntrial number of iterations for closed zero search - n number of dimension 4 or 6 - x intial guess for the closed orbit - tolx tolerance over the solution x - tolf tolerance over the evalution f(x) - - Output: - x closed orbit - - Return: - none - - Global variables: - status - - specific functions: - computeFandJS - ludcmp,lubksb - - Comments: - none - -****************************************************************************/ - -void Newton_RaphsonS(int ntrial, double x[], int n, double tolx) -{ - int k, i, *indx; - double errx, d, *bet, *fvect, **alpha; - - errx = 0.0; - // NR arrays start from 1 and not 0 !!! - indx = ivector(1, n); - bet = dvector(1, n); - fvect = dvector(1, n); - alpha = dmatrix(1, n, 1, n); - - for (k = 1; k <= ntrial; k++) { // loop over number of iterations - // supply function values at x in fvect and Jacobian matrix in fjac - computeFandJS(x, n, alpha, fvect); - - // Jacobian -Id - for (i = 1; i <= n; i++) - alpha[i][i] -= 1.0; - for (i = 1; i <= n; i++) - bet[i] = x[i] - fvect[i]; // right side of linear equation - // solve linear equations using LU decomposition using NR routines - dludcmp(alpha, n, indx, &d); - dlubksb(alpha, n, indx, bet); - errx = 0.0; // check root convergence - for (i = 1; i <= n; i++) { // update solution - errx += fabs(bet[i]); - x[i] += bet[i]; - } - - if (trace) - fprintf(stdout, - "%02d: cod % .5e % .5e % .5e % .5e % .5e % .5e errx =% .5e\n", - k, x[1], x[2], x[3], x[4], x[5], x[6], errx); - if (errx <= tolx) { - status.codflag = true; - break; - } - } - // check whever closed orbit found out - if ((k >= ntrial) && (errx >= tolx * 100)) status.codflag = false; - - free_dmatrix(alpha,1,n,1,n); free_dvector(bet,1,n); free_dvector(fvect,1,n); - free_ivector(indx,1,n); -} - - -/****************************************************************************/ -/* int Newton_Raphson(int n, double x[], int ntrial, double tolx) - - Purpose: - Newton_Rapson algorithm from Numerical Recipes - double precision algorithm - Robustess: quadratic convergence - Hint: for n-dimensional problem, the algo can be stuck on local minimum - In this case, it should be enough to provide a resonable starting - point. - - Method: - look for closed orbit solution of f(x) = x - This problems is equivalent to finding the zero of g(x)= f(x) - x - g(x+h) ~= f(x) - x + (Jacobian(f) -Id) h + O(h*h) - Then at first order we solve h: - h = - inverse(Jacobian(f) -Id) * (f(x)-x) - the new guess is then xnew = x + h - By iteration, this converges quadratically. - - The algo is stopped whenever |x -xnew| < tolx - - f(x) is computes by tracking over one turn - Jacobian(f) is computed numerically by numerical differentiation - These two operations are provided by the function computeFandJ - - Input: - ntrial number of iterations for closed zero search - x intial guess for the closed orbit - tolx tolerance over the solution x - tolf tolerance over the evalution f(x) - - Output: - x closed orbit - - Return: - none - - Global variables: - status - - specific functions: - computeFandJ - InvMat, LinTrans - - Comments: - none - -****************************************************************************/ -int Newton_Raphson (int n, Vector &x, int ntrial, double tolx) -{ - int k, i; - double errx; - Vector bet, fvect; - Matrix alpha; - - errx = 0.0; - - for (k = 1; k <= ntrial; k++) { // loop over number of iterations - // supply function values at x in fvect and Jacobian matrix in fjac - computeFandJ(n, x, alpha, fvect); - - // Jacobian - Id - for (i = 0; i < n; i++) - alpha[i][i] -= 1.0; - for (i = 0; i < n; i++) - bet[i] = x[i] - fvect[i]; // right side of linear equation - // inverse matrix using gauss jordan method from Tracy (from NR) - if (!InvMat((long) n,alpha)) - fprintf(stdout,"Matrix non inversible ...\n"); - LinTrans((long) n, alpha, bet); // bet = alpha*bet - errx = 0.0; // check root convergence - for (i = 0; i < n; i++) - { // update solution - errx += fabs(bet[i]); - x[i] += bet[i]; - } - - if (trace) - fprintf(stdout, - "%02d: cod2 % .5e % .5e % .5e % .5e % .5e % .5e errx =% .5e\n", - k, x[0], x[1], x[2], x[3], x[4], x[5], errx); - if (errx <= tolx) - { - status.codflag = true; - return 1; - } - } - // check whever closed orbit found out - if ((k >= ntrial) && (errx >= tolx)) - { - status.codflag = false; - return 1; - } - return 0; -} - - - -void rm_mean(long int n, double x[]) -{ - long int i; - double mean; - - mean = 0.0; - for (i = 0; i < n; i++) - mean += x[i]; - mean /= n; - for (i = 0; i < n; i++) - x[i] -= mean; -} +/* Tracy-2 + + J. Bengtsson, CBP, LBL 1990 - 1994 Pascal version + SLS, PSI 1995 - 1997 + M. Boege SLS, PSI 1998 C translation + L. Nadolski SOLEIL 2002 Link to NAFF, Radia field maps + J. Bengtsson NSLS-II, BNL 2004 - + +*/ + + +/**************************/ +/* Routines for printing */ +/**************************/ + +/**** same as asctime in C without the \n at the end****/ +char *asctime2(const struct tm *timeptr) +{ + // terminated with \0. + static char wday_name[7][4] = { + "Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat" + }; + // terminated with \0. + static char mon_name[12][4] = { + "Jan", "Feb", "Mar", "Apr", "May", "Jun", + "Jul", "Aug", "Sep", "Oct", "Nov", "Dec" + }; + static char result[26]; + + sprintf(result, "%.3s %.3s%3d %.2d:%.2d:%.2d %d", + wday_name[timeptr->tm_wday], + mon_name[timeptr->tm_mon], + timeptr->tm_mday, timeptr->tm_hour, + timeptr->tm_min, timeptr->tm_sec, + 1900 + timeptr->tm_year); + return result; +} + +/** Get time and date **/ +struct tm* GetTime() +{ + struct tm *whattime; + /* Get time and date */ + time_t aclock; + time(&aclock); /* Get time in seconds */ + whattime = localtime(&aclock); /* Convert time to struct */ + return whattime; +} + +/****************************************************************************/ +/* void printglob(void) + + Purpose: + Print global variables on screen + Print tunes and chromaticities + Print Oneturn matrix + + Input: + none + + Output: + output on the screen + + Return: + none + + Global variables: + globval + + Specific functions: + none + + Comments: + 26/03/03 Oneturn matrix added + 26/03/03 RF acceptance added + 10/05/03 Momentum compaction factor added + 16/05/03 Correction for a asymmetrical vaccum vessel + 20/06/03 Add corrector, skew quad and bpm number + 27/10/03 Add flag for radiation and chambre + + Comments copied from Tracy 2.7(soleil),Written by L.Nadolski. + +****************************************************************************/ +void printglob(void) +{ + printf("\n***************************************************************" + "***************\n"); + printf("\n"); + printf(" dPcommon = %9.3e dPparticle = %9.3e" + " Energy [GeV] = %.3f\n", + globval.dPcommon, globval.dPparticle, globval.Energy); + printf(" MaxAmplx [m] = %9.3e MaxAmply [m] = %9.3e" + " RFAccept [%%] = \xB1%4.2f\n", + Cell[0].maxampl[X_][0], Cell[0].maxampl[Y_][0], + globval.delta_RF*1e2); + printf(" MatMeth = %s ", globval.MatMeth ? "TRUE " : "FALSE"); + printf(" Cavity_On = %s ", globval.Cavity_on ? "TRUE " : "FALSE"); + printf(" Radiation_On = %s \n", + globval.radiation ? "TRUE " : "FALSE"); + printf(" bpm = %3ld qt = %3ld ", + globval.bpm, globval.qt); + printf(" Chambre_On = %s \n", status.chambre ? "TRUE " : "FALSE"); + printf(" hcorr = %3ld vcorr = %3ld\n\n", + globval.hcorr, globval.vcorr); + printf(" alphac = %8.4e\n", globval.Alphac); + printf(" nux = %13.6f nuz = %13.6f", + globval.TotalTune[X_], globval.TotalTune[Y_]); + if (globval.Cavity_on) + printf(" omega = %13.9f\n", globval.Omega); + else { + printf("\n"); + printf(" ksix = %13.6f ksiz = %13.6f\n", + globval.Chrom[X_], globval.Chrom[Y_]); + } + printf("\n"); + printf(" OneTurn matrix:\n"); + printf("\n"); + prtmat(ss_dim, globval.OneTurnMat); + fflush(stdout); +} + + + +/****************************************************************************/ +/* void printlatt(void) + + Purpose: + Print twiss parameters into file linlat.out + name, position, alpha, beta, phase, dispersion and its derivative + + Input: + none + + Output: + none + + Return: + none + + Global variables: + globval + + Specific functions: + getelem + + Comments: + 28/04/03 outfile end with .out extension instead of .dat + Twiss parameters are computed at the end of elements + +****************************************************************************/ +//void printlatt(void) +void printlatt(void) +{ + long int i = 0; + FILE *outf; + const char fic[] = "linlat.out"; + struct tm *newtime; + + /* Get time and date */ + newtime = GetTime(); + + if ((outf = fopen(fic,"w")) == NULL) + { + fprintf(stdout, "printlatt: Error while opening file %s \n",fic); + exit(1); + } + + fprintf(outf,"# TRACY III v. SYNCHROTRON SOLEIL -- %s -- %s \n", fic, asctime2(newtime)); + fprintf(outf, + "# name s alphax betax nux etax etapx alphay betay nuy etay etapy\n"); + fprintf(outf, + "# [m] [m] [m] [m] [m]\n"); + fprintf(outf, + "# exit \n"); + + for (i = 1L; i <= globval.Cell_nLoc; i++) + { + fprintf(outf, "%4ld:%.*s% 8.4f% 8.3f% 7.3f% 7.3f% 7.3f% 7.3f% 8.3f% 7.3f% 7.3f% 12.3e% 12.3e\n", + i, SymbolLength, Cell[i].Elem.PName, + Cell[i].S, Cell[i].Alpha[X_], Cell[i].Beta[X_], Cell[i].Nu[X_], Cell[i].Eta[X_], + Cell[i].Etap[X_], Cell[i].Alpha[Y_], Cell[i].Beta[Y_], Cell[i].Nu[Y_], + Cell[i].Eta[Y_], Cell[i].Etap[Y_]); + } + fclose(outf); +} + + + + + + +double int_curly_H1(long int n) +{ + /* Integration with Simpson's Rule */ + + double curly_H; + Vector2 alpha[3], beta[3], nu[3], eta[3], etap[3]; + + // only works for matrix style calculations + get_twiss3(n, alpha, beta, nu, eta, etap); + + curly_H = (get_curly_H(alpha[0][X_], beta[0][X_], eta[0][X_], etap[0][X_]) + +4.0*get_curly_H(alpha[1][X_], beta[1][X_], + eta[1][X_], etap[1][X_]) + + get_curly_H(alpha[2][X_], beta[2][X_], eta[2][X_], etap[2][X_])) + /6.0; + + return curly_H; +} + + +void prt_sigma(void) +{ + long int i; + double code = 0.0; + FILE *outf; + + outf = file_write("../out/sigma.out"); + + fprintf(outf, "# name s sqrt(sx) sqrt(sx') sqrt(sy) sqrt(sy')\n"); + fprintf(outf, "# [m] [mm] [mrad] [mm] [mrad]\n"); + fprintf(outf, "#\n"); + + for (i = 0; i <= globval.Cell_nLoc; i++) { + switch (Cell[i].Elem.Pkind) { + case drift: + code = 0.0; + break; + case Mpole: + if (Cell[i].Elem.M->Pirho != 0) + code = 0.5; + else if (Cell[i].Elem.M->PBpar[Quad+HOMmax] != 0) + code = sgn(Cell[i].Elem.M->PBpar[Quad+HOMmax]); + else if (Cell[i].Elem.M->PBpar[Sext+HOMmax] != 0) + code = 1.5*sgn(Cell[i].Elem.M->PBpar[Sext+HOMmax]); + else if (Cell[i].Fnum == globval.bpm) + code = 2.0; + else + code = 0.0; + break; + default: + code = 0.0; + break; + } + fprintf(outf, "%4ld %.*s %6.2f %4.1f %9.3e %9.3e %9.3e %9.3e\n", + i, SymbolLength, Cell[i].Elem.PName, Cell[i].S, code, + 1e3*sqrt(Cell[i].sigma[x_][x_]), + 1e3*sqrt(fabs(Cell[i].sigma[x_][px_])), + 1e3*sqrt(Cell[i].sigma[y_][y_]), + 1e3*sqrt(fabs(Cell[i].sigma[y_][py_]))); + } + + fclose(outf); +} + + +void recalc_S(void) +{ + long int k; + double S_tot; + + S_tot = 0.0; + for (k = 0; k <= globval.Cell_nLoc; k++) { + S_tot += Cell[k].Elem.PL; Cell[k].S = S_tot; + } +} + + +bool getcod(double dP, long &lastpos) +{ + return GetCOD(globval.CODimax, globval.CODeps, dP, lastpos); +} + + +void get_twiss3(long int loc, + Vector2 alpha[], Vector2 beta[], Vector2 nu[], + Vector2 eta[], Vector2 etap[]) +{ + /* Get Twiss functions at magnet entrance-, center-, and exit. */ + + long int j, k; + Vector2 dnu; + Matrix Ascr0, Ascr1; + + // Lattice functions at the magnet entrance + for (k = 0; k <= 1; k++) { + alpha[0][k] = Cell[loc-1].Alpha[k]; beta[0][k] = Cell[loc-1].Beta[k]; + nu[0][k] = Cell[loc-1].Nu[k]; + eta[0][k] = Cell[loc-1].Eta[k]; etap[0][k] = Cell[loc-1].Etap[k]; + } + + UnitMat(5, Ascr0); + for (k = 0; k <= 1; k++) { + Ascr0[2*k][2*k] = sqrt(Cell[loc-1].Beta[k]); Ascr0[2*k][2*k+1] = 0.0; + Ascr0[2*k+1][2*k] = -Cell[loc-1].Alpha[k]/Ascr0[2*k][2*k]; + Ascr0[2*k+1][2*k+1] = 1/Ascr0[2*k][2*k]; + } + Ascr0[0][4] = eta[0][X_]; Ascr0[1][4] = etap[0][X_]; + Ascr0[2][4] = eta[1][Y_]; Ascr0[3][4] = etap[1][Y_]; + CopyMat(5, Ascr0, Ascr1); MulLMat(5, Cell[loc].Elem.M->AU55, Ascr1); + dnu[0] = (atan(Ascr1[0][1]/Ascr1[0][0])-atan(Ascr0[0][1]/Ascr0[0][0])) + /(2.0*M_PI); + dnu[1] = (atan(Ascr1[2][3]/Ascr1[2][2])-atan(Ascr0[2][3]/Ascr0[2][2])) + /(2.0*M_PI); + + // Lattice functions at the magnet center + for (k = 0; k <= 1; k++) { + alpha[1][k] = -Ascr1[2*k][2*k]*Ascr1[2*k+1][2*k] + - Ascr1[2*k][2*k+1]*Ascr1[2*k+1][2*k+1]; + beta[1][k] = pow(Ascr1[2*k][2*k], 2.0) + pow(Ascr1[2*k][2*k+1], 2); + nu[1][k] = nu[0][k] + dnu[k]; + j = 2*k + 1; eta[1][k] = Ascr1[j-1][4]; etap[1][k] = Ascr1[j][4]; + } + + CopyMat(5, Ascr1, Ascr0); MulLMat(5, Cell[loc].Elem.M->AD55, Ascr1); + dnu[0] = (atan(Ascr1[0][1]/Ascr1[0][0])-atan(Ascr0[0][1]/Ascr0[0][0])) + /(2.0*M_PI); + dnu[1] = (atan(Ascr1[2][3]/Ascr1[2][2])-atan(Ascr0[2][3]/Ascr0[2][2])) + /(2.0*M_PI); + + // Lattice functions at the magnet exit + for (k = 0; k <= 1; k++) { + alpha[2][k] = -Ascr1[2*k][2*k]*Ascr1[2*k+1][2*k] + - Ascr1[2*k][2*k+1]*Ascr1[2*k+1][2*k+1]; + beta[2][k] = pow(Ascr1[2*k][2*k], 2.0) + pow(Ascr1[2*k][2*k+1], 2); + nu[2][k] = nu[1][k] + dnu[k]; + j = 2*k + 1; eta[2][k] = Ascr1[j-1][4]; etap[2][k] = Ascr1[j][4]; + } +} + + +void getabn(Vector2 &alpha, Vector2 &beta, Vector2 &nu) +{ + Vector2 gamma; + Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu); +} + + +void TraceABN(long i0, long i1, const Vector2 &alpha, const Vector2 &beta, + const Vector2 &eta, const Vector2 &etap, const double dP) +{ + long i, j; + double sb; + ss_vect<tps> Ascr; + + UnitMat(6, globval.Ascr); + for (i = 1; i <= 2; i++) { + sb = sqrt(beta[i-1]); j = i*2 - 1; + globval.Ascr[j-1][j-1] = sb; globval.Ascr[j-1][j] = 0.0; + globval.Ascr[j][j - 1] = -(alpha[i-1]/sb); globval.Ascr[j][j] = 1/sb; + } + globval.Ascr[0][4] = eta[0]; globval.Ascr[1][4] = etap[0]; + globval.Ascr[2][4] = eta[1]; globval.Ascr[3][4] = etap[1]; + + for (i = 0; i < 6; i++) + globval.CODvect[i] = 0.0; + globval.CODvect[4] = dP; + + if (globval.MatMeth) + Cell_Twiss_M(i0, i1, globval.Ascr, false, false, dP); + else { + for (i = 0; i <= 5; i++) { + Ascr[i] = tps(globval.CODvect[i]); + for (j = 0; j <= 5; j++) + Ascr[i] += globval.Ascr[i][j]*tps(0.0, j+1); + } + Cell_Twiss(i0, i1, Ascr, false, false, dP); + } + +} + + +void FitTune(long qf, long qd, double nux, double nuy) +{ + long i; + iVector2 nq = {0,0}; + Vector2 nu = {0.0, 0.0}; + fitvect qfbuf, qdbuf; + + /* Get elements for the first quadrupole family */ + nq[X_] = GetnKid(qf); + for (i = 1; i <= nq[X_]; i++) + qfbuf[i-1] = Elem_GetPos(qf, i); + + /* Get elements for the second quadrupole family */ + nq[Y_] = GetnKid(qd); + for (i = 1; i <= nq[Y_]; i++) + qdbuf[i - 1] = Elem_GetPos(qd, i); + + nu[X_] = nux; nu[Y_] = nuy; + + /* fit tunes */ + Ring_Fittune(nu, nueps, nq, qfbuf, qdbuf, nudkL, nuimax); +} + + +void FitChrom(long sf, long sd, double ksix, double ksiy) +{ + long i; + iVector2 ns = {0,0}; + fitvect sfbuf, sdbuf; + Vector2 ksi ={0.0, 0.0}; + + /* Get elements for the first sextupole family */ + ns[X_] = GetnKid(sf); + for (i = 1; i <= ns[X_]; i++) + sfbuf[i-1] = Elem_GetPos(sf, i); + + /* Get elements for the second sextupole family */ + ns[Y_] = GetnKid(sd); + for (i = 1; i <= ns[Y_]; i++) + sdbuf[i-1] = Elem_GetPos(sd, i); + + ksi[X_] = ksix; ksi[Y_] = ksiy; + + /* Fit chromaticities */ + /* Ring_Fitchrom(ksi, ksieps, ns, sfbuf, sdbuf, 1.0, 1);*/ + Ring_Fitchrom(ksi, ksieps, ns, sfbuf, sdbuf, ksidkpL, ksiimax); +} + + +void FitDisp(long q, long pos, double eta) +{ + long i, nq; + fitvect qbuf; + + /* Get elements for the quadrupole family */ + nq = GetnKid(q); + for (i = 1; i <= nq; i++) + qbuf[i-1] = Elem_GetPos(q, i); + + Ring_FitDisp(pos, eta, dispeps, nq, qbuf, dispdkL, dispimax); +} + + +#define nfloq 4 + +void getfloqs(Vector &x) +{ + // Transform to Floquet space + LinTrans(nfloq, globval.Ascrinv, x); +} + +#undef nfloq + + +#define ntrack 4 + +// 4D tracking in normal or Floquet space over nmax turns + +void track(const char *file_name, + double ic1, double ic2, double ic3, double ic4, double dp, + long int nmax, long int &lastn, long int &lastpos, int floqs, + double f_rf) +{ + /* Single particle tracking around closed orbit: + + Output floqs + + Phase Space 0 [x, px, y, py, delta, ct] + Floquet Space 1 [x^, px^, y^, py^, delta, ct] + Action-Angle Variables 2 [2Jx, phx, 2Jy, phiy, delta, ct] + + */ + long int i; + double twoJx, twoJy, phix, phiy, scl_1 = 1.0, scl_2 = 1.0; + Vector x0, x1, x2, xf; + FILE *outf; + + bool prt = false; + + if ((floqs == 0)) { + scl_1 = 1e3; scl_2 = 1e3; + x0[x_] = ic1; x0[px_] = ic2; x0[y_] = ic3; x0[py_] = ic4; + } else if ((floqs == 1)) { + scl_1 = 1.0; scl_2 = 1.0; + x0[x_] = ic1; x0[px_] = ic2; x0[y_] = ic3; x0[py_] = ic4; + LinTrans(4, globval.Ascr, x0); + } else if (floqs == 2) { + scl_1 = 1e6; scl_2 = 1.0; + x0[x_] = sqrt(ic1)*cos(ic2); x0[px_] = -sqrt(ic1)*sin(ic2); + x0[y_] = sqrt(ic3)*cos(ic4); x0[py_] = -sqrt(ic3)*sin(ic4); + LinTrans(4, globval.Ascr, x0); + } + + outf = file_write(file_name); + fprintf(outf, "# Tracking with TRACY"); + getcod(dp, lastpos); + if (floqs == 0) + fprintf(outf, "\n"); + else if (floqs == 1) { + Ring_GetTwiss(false, dp); + fprintf(outf, "# (Floquet space)\n"); + } else if (floqs == 2) { + Ring_GetTwiss(false, dp); + fprintf(outf, "# (Action-Angle variables)\n"); + } + fprintf(outf, "#\n"); + fprintf(outf, "#%3d%6ld% .1E% .1E% .1E% .1E% 7.5f% 7.5f\n", + 1, nmax, 1e0, 1e0, 0e0, 0e0, globval.TotalTune[0], + globval.TotalTune[1]); + if (floqs == 0) { + fprintf(outf, "# N x p_x y p_y"); + fprintf(outf, " delta cdt\n"); + fprintf(outf, "# [mm] [mrad]" + " [mm] [mrad]"); + } else if (floqs == 1) { + fprintf(outf, "# N x^ px^ y^ py^"); + fprintf(outf, " delta cdt"); + fprintf(outf, "# " + " "); + } else if (floqs == 2) { + fprintf(outf, "# N 2Jx phi_x 2Jy phi_y"); + fprintf(outf, " delta cdt\n"); + fprintf(outf, "# " + " "); + } + if (f_rf == 0.0){ + fprintf(outf, " [%%] [mm]\n"); + fprintf(outf, "#\n"); + fprintf(outf, "%4d %23.16e %23.16e %23.16e %23.16e %23.16e %23.16e\n", + 0, scl_1*ic1, scl_2*ic2, scl_1*ic3, scl_2*ic4, 1e2*dp, + 1e3*globval.CODvect[ct_]); + } else { + fprintf(outf, " [%%] [deg]\n"); + fprintf(outf, "#\n"); + fprintf(outf, "%4d %23.16e %23.16e %23.16e %23.16e %23.16e %23.16e\n", + 0, scl_1*ic1, scl_2*ic2, scl_1*ic3, scl_2*ic4, 1e2*dp, + 2.0*f_rf*180.0*globval.CODvect[ct_]/c0); + } + x2[x_] = x0[x_] + globval.CODvect[x_]; + x2[px_] = x0[px_] + globval.CODvect[px_]; + x2[y_] = x0[y_] + globval.CODvect[y_]; + x2[py_] = x0[py_] + globval.CODvect[py_]; + if (globval.Cavity_on) { + x2[delta_] = dp + globval.CODvect[delta_]; x2[ct_] = globval.CODvect[ct_]; + } else { + x2[delta_] = dp; x2[ct_] = 0.0; + } + + lastn = 0; + + if (prt) { + printf("\n"); + printf("track:\n"); + printf("%4ld %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", + lastn, 1e3*x2[x_], 1e3*x2[px_], 1e3*x2[y_], 1e3*x2[py_], + 1e2*x2[delta_], 1e3*x2[ct_]); + } + + if (globval.MatMeth) Cell_Concat(dp); + + do { + (lastn)++; + for (i = 0; i < nv_; i++) + x1[i] = x2[i]; + + if (globval.MatMeth) { + Cell_fPass(x2, lastpos); + } else { + Cell_Pass(0, globval.Cell_nLoc, x2, lastpos); + } + + for (i = x_; i <= py_; i++) + xf[i] = x2[i] - globval.CODvect[i]; + + for (i = delta_; i <= ct_; i++) + if (globval.Cavity_on && (i != ct_)) + xf[i] = x2[i] - globval.CODvect[i]; + else + xf[i] = x2[i]; + + if (floqs == 1) + getfloqs(xf); + else if (floqs == 2) { + getfloqs(xf); + twoJx = pow(xf[x_], 2) + pow(xf[px_], 2); + twoJy = pow(xf[y_], 2) + pow(xf[py_], 2); + phix = atan2(xf[px_], xf[x_]); + phiy = atan2(xf[py_], xf[y_]); + xf[x_] = twoJx; xf[px_] = phix; xf[y_] = twoJy; xf[py_] = phiy; + } + if (f_rf == 0.0) + fprintf(outf, + "%4ld %23.16le %23.16le %23.16le %23.16le %23.16le %23.16le\n", + lastn, scl_1*xf[0], scl_2*xf[1], scl_1*xf[2], scl_2*xf[3], + 1e2*xf[4], 1e3*xf[5]); + else + fprintf(outf, + "%4ld %23.16le %23.16le %23.16le %23.16le %23.16le %23.16le\n", + lastn, scl_1*xf[0], scl_2*xf[1], scl_1*xf[2], scl_2*xf[3], + 1e2*xf[4], 2.0*f_rf*180.0*xf[5]/c0); + } while ((lastn != nmax) && (lastpos == globval.Cell_nLoc)); + + if (globval.MatMeth) + Cell_Pass(0, globval.Cell_nLoc, x1, lastpos); + + fclose(outf); +} + +#undef ntrack + + +#define step 0.1 +#define px 0.0 +#define py 0.0 +void track_(double r, struct LOC_getdynap *LINK) +{ + long i, lastn, lastpos; + Vector x; + + x[0] = r * cos(LINK->phi); + x[1] = px; + x[2] = r * sin(LINK->phi); + x[3] = py; + x[4] = LINK->delta; + x[5] = 0.0; + /* transform to phase space */ + if (LINK->floqs) { + LinTrans(5, globval.Ascr, x); + } + for (i = 0; i <= 3; i++) + x[i] += globval.CODvect[i]; + lastn = 0; + do { + lastn++; + if (globval.MatMeth) { + Cell_fPass(x, lastpos); + } else { + Cell_Pass(0, globval.Cell_nLoc, x, lastpos); + } + } while (lastn != LINK->nturn && lastpos == globval.Cell_nLoc); + LINK->lost = (lastn != LINK->nturn); +} +#undef step +#undef px +#undef py + + +/****************************************************************************/ +/* void Trac(double x, double px, double y, double py, double dp, double ctau, + long nmax, long pos, long &lastn, long &lastpos, FILE *outf1) + + Purpose: + Single particle tracking w/ respect to the chamber centrum + + Input: + x, px, y, py 4 transverses coordinates + ctau c*tau + dp energy offset + nmax number of turns + pos starting position for tracking + aperture global physical aperture + + Output: + lastn last n (should be nmax if not lost) + lastpos last position in the ring + outf1 output file with 6D phase at different pos + + Return: + none + + Global variables: + globval + + specific functions: + Cell_Pass + + Comments: + Absolute TRACKING with respect to the center of the vacuum vessel + BUG: last printout is wrong because not at pos but at the end of + the ring + 26/04/03 print output for phase space is for position pos now + 01/12/03 tracking from 0 to pos -1L instead of 0 to pos + (wrong observation point) + +****************************************************************************/ +void Trac(double x, double px, double y, double py, double dp, double ctau, + long nmax, long pos, long &lastn, long &lastpos, FILE *outf1) +{ + bool lostF; /* Lost particle Flag */ + Vector x1; /* tracking coordinates */ + Vector2 aperture; + + /* Compute closed orbit : usefull if insertion devices */ + + aperture[0] = 1e0; + aperture[1] = 1e0; + + x1[0] = x; x1[1] = px; + x1[2] = y; x1[3] = py; + x1[4] =dp; x1[5] = ctau; + + lastn = 0L; + lostF = true; + + (lastpos)=pos; + if(trace) fprintf(outf1, "\n"); + fprintf(outf1, "%6ld %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e \n", lastn, + x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]); + Cell_Pass(pos -1L, globval.Cell_nLoc, x1, lastpos); + + if (lastpos == globval.Cell_nLoc) + do { + (lastn)++; + Cell_Pass(0L,pos-1L, x1, lastpos); + if(!trace) { + fprintf(outf1, "%6ld %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e \n", + lastn, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]); + } + if (lastpos == pos-1L) Cell_Pass(pos-1L,globval.Cell_nLoc, x1, lastpos); + } + while (((lastn) < nmax) && ((lastpos) == globval.Cell_nLoc)); + + if (lastpos == globval.Cell_nLoc) Cell_Pass(0L,pos, x1, lastpos); + + if (lastpos != pos) { + printf("Trac: Particle lost \n"); + fprintf(stdout, "turn:%6ld plane: %1d" + " %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", + lastn, status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]); + } +} + +/****************************************************************************/ +/*bool chk_if_lost(double x0, double y0, double delta, + long int nturn, bool floqs) + + Purpose: + Binary search for dynamical aperture in Floquet space. + + Input: + none + + Output: + none + + Return: + none + + Global variables: + px_0, py_0 + + Specific functions: + chk_if_lost + + Comments: + none + +****************************************************************************/ + +#define nfloq 4 +bool chk_if_lost(double x0, double y0, double delta, + long int nturn, bool floqs) +{ + long int i, lastn, lastpos; + Vector x; + + bool prt = false; + + x[x_] = x0; x[px_] = px_0; x[y_] = y0; x[py_] = py_0; + x[delta_] = delta; x[ct_] = 0.0; + if (floqs) + // transform to phase space + LinTrans(nfloq, globval.Ascr, x); + for (i = 0; i <= 3; i++) + x[i] += globval.CODvect[i]; + + lastn = 0; + if (prt) { + printf("\n"); + printf("chk_if_lost:\n"); + printf("%4ld %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", + lastn, 1e3*x[x_], 1e3*x[px_], 1e3*x[y_], 1e3*x[py_], + 1e2*x[delta_], 1e3*x[ct_]); + } + do { + lastn++; + if (globval.MatMeth) + Cell_fPass(x, lastpos); + else + Cell_Pass(0, globval.Cell_nLoc, x, lastpos); + if (prt) + printf("%4ld %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", + lastn, 1e3*x[x_], 1e3*x[px_], 1e3*x[y_], 1e3*x[py_], + 1e2*x[delta_], 1e3*x[ct_]); + } while ((lastn != nturn) && (lastpos == globval.Cell_nLoc)); + return(lastn != nturn); +} +#undef nfloq + +/****************************************************************************/ +/* void getdynap(double *r, double phi, double delta, double eps, + int nturn, bool floqs) + + Purpose: + Binary search for dynamical aperture in Floquet space. + + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + chk_if_lost + + Comments: + none + +****************************************************************************/ +void getdynap(double &r, double phi, double delta, double eps, + int nturn, bool floqs) +{ + /* Determine dynamical aperture by binary search. */ + double rmin = 0.0, rmax = r; + + const bool prt = false; + const double r_reset = 1e-3, r0 = 10e-3; + + if (prt) printf("\n"); + + if (globval.MatMeth) Cell_Concat(delta); + while (!chk_if_lost(rmax*cos(phi), rmax*sin(phi), delta, nturn, floqs)) { + if (rmax < r_reset) rmax = r0; + rmax *= 2.0; + } + while (rmax-rmin >= eps) { + r = rmin + (rmax-rmin)/2.0; + if (prt) printf("getdynap: %6.3f %6.3f %6.3f\n", + 1e3*rmin, 1e3*rmax, 1e3*r); + if (! chk_if_lost(r*cos(phi), r*sin(phi), delta, nturn, floqs) ) + rmin = r; + else + rmax = r; + if (rmin > rmax) { + printf("getdynap: rmin > rmax\n"); + exit_(0); + } + + } + r = rmin; +} + + + +/****************************************************************************/ +/* void getcsAscr(void) + + Purpose: + Get Courant-Snyder Ascr + + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + none + +****************************************************************************/ +void getcsAscr(void) +{ + long i, j; + double phi; + Matrix R; + + UnitMat(6, R); + for (i = 1; i <= 2; i++) { + phi = -atan2(globval.Ascr[i * 2 - 2][i * 2 - 1], globval.Ascr[i * 2 - 2] + [i * 2 - 2]); + R[i * 2 - 2][i * 2 - 2] = cos(phi); + R[i * 2 - 1][i * 2 - 1] = R[i * 2 - 2][i * 2 - 2]; + R[i * 2 - 2][i * 2 - 1] = sin(phi); + R[i * 2 - 1][i * 2 - 2] = -R[i * 2 - 2][i * 2 - 1]; + } + MulRMat(6, globval.Ascr, R); + for (i = 1; i <= 2; i++) { + if (globval.Ascr[i * 2 - 2][i * 2 - 2] < 0.0) { + for (j = 0; j <= 5; j++) { + globval.Ascr[j][i * 2 - 2] = -globval.Ascr[j][i * 2 - 2]; + globval.Ascr[j][i * 2 - 1] = -globval.Ascr[j][i * 2 - 1]; + } + } + } + if (!InvMat(6, globval.Ascrinv)) + printf(" *** Ascr is singular\n"); +} + + +/****************************************************************************/ +/* void dynap(double r, double delta, double eps, int npoint, int nturn, + double x[], double y[], bool floqs, bool print) + + Purpose: + Determine the dynamical aperture by tracking using polar coordinates, + and sampling in phase. + Assumes mid-plane symmetry + + Input: + r initial guess + delta off momentum energy + eps precision for binary search + npoint sample number for phase coordinate + nturn number of turn for computing da + floqs true means Floquet space + print true means Print out to the screen + + Output: + x[] horizontal dynamics aperture + y[] vertical dynamics aperture + + Return: + none + + Global variables: + none + + Specific functions: + getdynap + + Comments: + none + +****************************************************************************/ +void dynap(FILE *fp, double r, const double delta, + const double eps, const int npoint, const int nturn, + double x[], double y[], const bool floqs, const bool print) + +{ + /* Determine the dynamical aperture by tracking. + Assumes mid-plane symmetry. */ + + long int i, lastpos; + double phi, x_min, x_max, y_min, y_max; + + getcod(delta, lastpos); + if (floqs) { + Ring_GetTwiss(false, delta); + if (print) { + printf("\n"); + printf("Dynamical Aperture (Floquet space):\n"); + printf(" x^ y^\n"); + printf("\n"); + } + fprintf(fp, "# Dynamical Aperture (Floquet space):\n"); + fprintf(fp, "# x^ y^\n"); + fprintf(fp, "#\n"); + } else { + if (print) { + printf("\n"); + printf("Dynamical Aperture:\n"); + printf(" x y\n"); + printf(" [mm] [mm]\n"); + printf("\n"); + } + fprintf(fp, "# Dynamical Aperture:\n"); + fprintf(fp, "# x y\n"); + fprintf(fp, "# [mm] [mm]\n"); + fprintf(fp, "#\n"); + } + + x_min = 0.0; x_max = 0.0; y_min = 0.0; y_max = 0.0; + for (i = 0; i < npoint; i++) { + phi = i*M_PI/(npoint-1); + if (i == 0) + phi = 1e-3; + else if (i == npoint-1) + phi -= 1e-3; + getdynap(r, phi, delta, eps, nturn, floqs); + x[i] = r*cos(phi); y[i] = r*sin(phi); + x_min = min(x[i], x_min); x_max = max(x[i], x_max); + y_min = min(y[i], y_min); y_max = max(y[i], y_max); + if (!floqs) { + if (print) + printf(" %6.2f %6.2f\n", 1e3*x[i], 1e3*y[i]); + fprintf(fp, " %6.2f %6.2f\n", 1e3*x[i], 1e3*y[i]); + } else { + if (print) + printf(" %10.3e %10.3e\n", x[i], y[i]); + fprintf(fp, " %10.3e %10.3e\n", x[i], y[i]); + } + fflush(fp); + } + + if (print) { + printf("\n"); + printf(" x^: %6.2f - %5.2f y^: %6.2f - %5.2f mm\n", + 1e3*x_min, 1e3*x_max, 1e3*y_min, 1e3*y_max); + } +} + +/****************************************************************************/ +/* double get_aper(int n, double x[], double y[]) + + Purpose: + + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + none + +****************************************************************************/ +double get_aper(int n, double x[], double y[]) +{ + int i; + double A; + + A = 0.0; + for (i = 2; i <= n; i++) + A += x[i-2]*y[i-1] - x[i-1]*y[i-2]; + A += x[n-1]*y[0] - x[0]*y[n-1]; +// x2 from mid-plane symmetry + A = fabs(A); +// printf("\n"); +// printf(" Dyn. Aper.: %5.1f mm^2\n", 1e6*A); + return(A); +} + + +void GetTrack(const char *file_name, + long *n, double x[], double px[], double y[], double py[]) +{ + int k; + char line[200]; + FILE *inf; + + inf = file_read(file_name); + + do { + fgets(line, 200, inf); + } while (strstr(line, "#") != NULL); + + // skip initial conditions + fgets(line, 200, inf); + + do { + sscanf(line, "%d", &k); + sscanf(line, "%*d %lf %lf %lf %lf", &x[k-1], &px[k-1], &y[k-1], &py[k-1]); + } while (fgets(line, 200, inf) != NULL); + + *n = k; + + fclose(inf); +} + + +/****************************************************************************/ +/* void Getj(long n, double *x, double *px, double *y, double *py) + + Purpose: + Calculates the linear invariant + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + none + +****************************************************************************/ +void Getj(long n, double *x, double *px, double *y, double *py) +{ + long int i; + + for (i = 0; i < n; i++) { + x[i] = (pow(x[i], 2)+pow(px[i], 2))/2.0; + y[i] = (pow(y[i], 2)+pow(py[i], 2))/2.0; + } +} + +/****************************************************************************/ +/* double GetArg(double x, double px, double nu) + + Purpose: + get argument of x + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + 17/07/03 use M_PI instead of pi + +****************************************************************************/ +double GetArg(double x, double px, double nu) +{ + double phi, val; + + phi = GetAngle(x, px); + if (phi < 0.0) + phi += 2.0 * M_PI; + val = phi + Fract(nu) * 2.0 * M_PI; + if (val < 0.0) + val += 2.0 * M_PI; + return val; +} + +/****************************************************************************/ +/* void GetPhi(long n, double *x, double *px, double *y, double *py) + + Purpose: + get linear phases + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + none + +****************************************************************************/ +void GetPhi(long n, double *x, double *px, double *y, double *py) +{ + /* Calculates the linear phase */ + long i; + + for (i = 1; i <= n; i++) { + x[i - 1] = GetArg(x[i - 1], px[i - 1], i * globval.TotalTune[0]); + y[i - 1] = GetArg(y[i - 1], py[i - 1], i * globval.TotalTune[1]); + } +} + + +/*********************************/ +/* Routines for Fourier analysis */ +/*********************************/ + +void Sinfft(int n, double xr[]) +{ + /* DFT with sine window */ + int i; + double xi[n]; + + for (i = 0; i < n; i++) { + xr[i] = sin((double)i/(double)n*M_PI)*xr[i]; xi[i] = 0.0; + } + FFT(n, xr, xi); + for (i = 0; i < n; i++) + xr[i] = sqrt(xr[i]*xr[i]+xi[i]*xi[i]); +} + +void sin_FFT(int n, double xr[]) +{ + /* DFT with sine window */ + int i; + double *xi; + + xi = dvector(1, 2*n); + + for (i = 1; i <= n; i++) { + xi[2*i-1] = sin((double)i/n*M_PI)*xr[i-1]; xi[2*i] = 0.0; + } + dfour1(xi, (unsigned long)n, 1); + for (i = 1; i <= n; i++) + xr[i-1] = sqrt(pow(xi[2*i-1], 2)+pow(xi[2*i], 2))*2.0/n; + + free_dvector(xi, 1, 2*n); +} + + +void sin_FFT(int n, double xr[], double xi[]) +{ + /* DFT with sine window */ + int i; + double *xri; + + xri = dvector(1, 2*n); + + for (i = 1; i <= n; i++) { + xri[2*i-1] = sin((double)i/n*M_PI)*xr[i-1]; + xri[2*i] = sin((double)i/n*M_PI)*xi[i-1]; + } + dfour1(xri, (unsigned long)n, 1); + for (i = 1; i <= n; i++) { + xr[i-1] = sqrt(pow(xri[2*i-1], 2)+pow(xri[2*i], 2))*2.0/n; + xi[i-1] = atan2(xri[2*i], xri[2*i-1]); + } + + free_dvector(xri, 1, 2*n); +} + + +void GetInd(int n, int k, int *ind1, int *ind3) +{ + if (k == 1) { + *ind1 = 2; *ind3 = 2; + } else if (k == n/2+1) { + *ind1 = n/2; *ind3 = n/2; + } else { + *ind1 = k - 1; *ind3 = k + 1; + } +} + + +void GetInd1(int n, int k, int *ind1, int *ind3) +{ + if (k == 1) { + *ind1 = 2; *ind3 = 2; + } else if (k == n) { + *ind1 = n - 1; *ind3 = n - 1; + } else { + *ind1 = k - 1; *ind3 = k + 1; + } +} + + +void GetPeak(int n, double *x, int *k) +{ + /* Locate peak in DFT spectrum */ + int ind1, ind2, ind3; + double peak; + + peak = 0.0; *k = 1; + for (ind2 = 1; ind2 <= n/2+1; ind2++) { + GetInd(n, ind2, &ind1, &ind3); + if (x[ind2-1] > peak && x[ind1-1] < x[ind2-1] && + x[ind3-1] < x[ind2-1]) + { + peak = x[ind2-1]; *k = ind2; + } + } +} + + +void GetPeak1(int n, double *x, int *k) +{ + /* Locate peak in DFT spectrum */ + int ind1, ind2, ind3; + double peak; + + peak = 0.0; *k = 1; + for (ind2 = 1; ind2 <= n; ind2++) { + GetInd1(n, ind2, &ind1, &ind3); + if (x[ind2-1] > peak && x[ind1-1] < x[ind2-1] && + x[ind3-1] < x[ind2-1]) + { + peak = x[ind2-1]; *k = ind2; + } + } +} + + +double Int2snu(int n, double *x, int k) +{ + /* Get frequency by nonlinear interpolation with two samples + for sine window. The interpolation is: + + 1 2 A(k) 1 + nu = - [ k - 1 + ------------- - - ] , k-1 <= N nu <= k + N A(k-1) + A(k) 2 + */ + int ind, ind1, ind3; + double ampl1, ampl2; + + GetInd(n, k, &ind1, &ind3); + if (x[ind3 - 1] > x[ind1 - 1]) { + ampl1 = x[k - 1]; + ampl2 = x[ind3 - 1]; + ind = k; + } else { + ampl1 = x[ind1 - 1]; + ampl2 = x[k - 1]; + /* Interpolate in right direction for 0 frequency */ + if (k != 1) + ind = ind1; + else + ind = 0; + } + /* Avoid division by zero */ + if (ampl1 + ampl2 != 0.0) + return ((ind - 1 + 2 * ampl2 / (ampl1 + ampl2) - 0.5) / n); + else + return 0.0; +} + + +/****************************************************************************/ +/* double Sinc(double omega) + + Purpose: + + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + zero test to be changed numericallywise + +****************************************************************************/ +double Sinc(double omega) +{ + /* Function to calculate: + + sin( omega ) + ------------ + omega + */ + if (omega != 0.0) + return (sin(omega) / omega); + else + return 1.0; +} + + +/****************************************************************************/ +/* double intsampl(long n, double *x, double nu, long k) + + Purpose: + + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + 17/07/03 use M_PI instead of pi + +****************************************************************************/ +double intsampl(int n, double *x, double nu, int k) +{ + /* Get amplitude by nonlinear interpolation for sine window. The + distribution is given by: + + 1 sin pi ( k + 1/2 ) sin pi ( k - 1/2 ) + F(k) = - ( -------------------- + -------------------- ) + 2 pi ( k + 1/2 ) pi ( k - 1/2 ) + */ + double corr; + + corr = (Sinc(M_PI * (k - 1 + 0.5 - nu * n)) + + Sinc(M_PI * (k - 1 - 0.5 - nu * n))) / 2; + return (x[k - 1] / corr); +} + + +/****************************************************************************/ +/* double linint(long n, long k, double nu, double *x) + + Purpose: + + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + 17/07/03 use M_PI instead of pi + +****************************************************************************/ +double linint(int n, int k, double nu, double *x) +{ + /* Get phase by linear interpolation for rectangular window + with -pi <= phi <= pi */ + int i; + double phi; + double xr[n], xi[n]; + + for (i = 0; i < n; i++) { + xr[i] = x[i]; xi[i] = 0.0; + } + FFT(n, xr, xi); + phi = GetAngle(xr[k-1], xi[k-1]) - (n*nu-k+1)*M_PI; + if (phi > M_PI) + phi -= 2.0*M_PI; + else if (phi < -M_PI) + phi += 2.0*M_PI; + return phi; +} + +/****************************************************************************/ +/* void FndRes(struct LOC_findres *LINK) + + Purpose: + + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + none + +****************************************************************************/ +void FndRes(struct LOC_findres *LINK) +{ + int i, j, FORLIM, FORLIM1; + double delta; + + FORLIM = LINK->n; + for (i = 0; i <= FORLIM; i++) { + FORLIM1 = LINK->n; + for (j = -LINK->n; j <= FORLIM1; j++) { + delta = fabs(i * LINK->nux + j * LINK->nuy); + delta -= (int)delta; + if (delta > 0.5) + delta = 1 - delta; + delta = fabs(delta - LINK->f); + delta -= (int)delta; + if (delta > 0.5) + delta = 1 - delta; + if (delta < LINK->eps) { + if (abs(i) + abs(j) < LINK->n && (i != 0 || j >= 0)) { + LINK->found = true; + *LINK->nx = i; + *LINK->ny = j; + } + } + } + } +} + + +/****************************************************************************/ +/* void FindRes(long n_, double nux_, double nuy_, double f_, + long *nx_, long *ny_) + + Purpose: + + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + none + +****************************************************************************/ +void FindRes(int n_, double nux_, double nuy_, double f_, int *nx_, int *ny_) +{ + /* Match f by a linear combination of nux and nuy */ + struct LOC_findres V; + + V.n = n_; + V.nux = nux_; + V.nuy = nuy_; + V.f = f_; + V.nx = nx_; + V.ny = ny_; + V.found = false; + V.eps = 0.5e-6; + do { + V.eps = 10 * V.eps; + FndRes(&V); + } while (!V.found); +} + + +void GetPeaks(int n, double *x, int nf, double *nu, double *A) +{ + int i, k, ind1, ind3; + + for (i = 0; i < nf; i++) { + GetPeak(n, x, &k); + nu[i] = Int2snu(n, x, k); A[i] = intsampl(n, x, nu[i], k); + /* Make peak flat to allow for new call */ + GetInd(n, k, &ind1, &ind3); + if (x[ind1-1] > x[ind3-1]) + x[k-1] = x[ind1-1]; + else + x[k-1] = x[ind3-1]; + } +} + + +void GetPeaks1(int n, double *x, int nf, double *nu, double *A) +{ + int i, k, ind1, ind3; + + for (i = 0; i < nf; i++) { + GetPeak1(n, x, &k); + nu[i] = Int2snu(n, x, k); A[i] = intsampl(n, x, nu[i], k); + /* Make peak flat to allow for new call */ + GetInd1(n, k, &ind1, &ind3); + if (x[ind1-1] > x[ind3-1]) + x[k-1] = x[ind1-1]; + else + x[k-1] = x[ind3-1]; + } +} + +/*******************************/ +/* Routines for magnetic error */ +/*******************************/ + +void SetTol(int Fnum, double dxrms, double dyrms, double drrms) +{ + int i; + long k; + + for (i = 1; i <= GetnKid(Fnum); i++) { + k = Elem_GetPos(Fnum, i); + Cell[k].Elem.M->PdSrms[X_] = dxrms; + Cell[k].Elem.M->PdSrnd[X_] = normranf(); + Cell[k].Elem.M->PdSrms[Y_] = dyrms; + Cell[k].Elem.M->PdSrnd[Y_] = normranf(); + Cell[k].Elem.M->PdTrms = drrms; + Cell[k].Elem.M->PdTrnd = normranf(); + Mpole_SetdS(Fnum, i); Mpole_SetdT(Fnum, i); + } +} + + +void Scale_Tol(int Fnum, double dxrms, double dyrms, double drrms) +{ + int Knum; + long int loc; + + for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) { + loc = Elem_GetPos(Fnum, Knum); + Cell[loc].Elem.M->PdSrms[X_] = dxrms; Cell[loc].Elem.M->PdSrms[Y_] = dyrms; + Cell[loc].Elem.M->PdTrms = drrms; + Mpole_SetdS(Fnum, Knum); Mpole_SetdT(Fnum, Knum); + } +} + + +/****************************************************************************/ +/* void SetaTol(int Fnum, int Knum, double dx, double dy, double dr) + + Purpose: + Set a known random multipole displacement error + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + none + +****************************************************************************/ +void SetaTol(int Fnum, int Knum, double dx, double dy, double dr) +{ + long int loc; + + loc = Elem_GetPos(Fnum, Knum); + Cell[loc].Elem.M->PdSrms[0] = dx; Cell[loc].Elem.M->PdSrnd[0] = 1e0; + Cell[loc].Elem.M->PdSrms[1] = dy; Cell[loc].Elem.M->PdSrnd[1] = 1e0; + Cell[loc].Elem.M->PdTrms = dr; Cell[loc].Elem.M->PdTrnd = 1e0; + Mpole_SetdS(Fnum, Knum); Mpole_SetdT(Fnum, Knum); +} + + +void ini_aper(const double Dxmin, const double Dxmax, + const double Dymin, const double Dymax) +{ + int k; + + for (k = 0; k <= globval.Cell_nLoc; k++) { + Cell[k].maxampl[X_][0] = Dxmin; Cell[k].maxampl[X_][1] = Dxmax; + Cell[k].maxampl[Y_][0] = Dymin; Cell[k].maxampl[Y_][1] = Dymax; + } +} + +void set_aper(const int Fnum, const double Dxmin, const double Dxmax, + const double Dymin, const double Dymax) +{ + int i; + long int loc; + + for (i = 1; i <= GetnKid(Fnum); i++) { + loc = Elem_GetPos(Fnum, i); + Cell[loc].maxampl[X_][0] = Dxmin; Cell[loc].maxampl[X_][1] = Dxmax; + Cell[loc].maxampl[Y_][0] = Dymin; Cell[loc].maxampl[Y_][1] = Dymax; + } +} + + +void LoadApertures(const char *ChamberFileName) +{ + char line[128], FamName[32]; + long Fnum; + double Xmin, Xmax, Ymin, Ymax; + FILE *ChamberFile; + + ChamberFile = file_read(ChamberFileName); + + do + fgets(line, 128, ChamberFile); + while (strstr(line, "#") != NULL); + + do { + sscanf(line,"%s %lf %lf %lf %lf", FamName,&Xmin, &Xmax, &Ymin,&Ymax); + Fnum = ElemIndex(FamName); + if (Fnum > 0) set_aper(Fnum, Xmin, Xmax, Ymin, Ymax); + } while (fgets(line, 128, ChamberFile ) != NULL); + + fclose(ChamberFile); +} + + +// Load tolerances from the file +void LoadTolerances(const char *TolFileName) +{ + char line[128], FamName[32]; + int Fnum; + double dx, dy, dr; + FILE *tolfile; + + tolfile = file_read(TolFileName); + + do + fgets(line, 128, tolfile); + while (strstr(line, "#") != NULL); + + do { + if (strstr(line, "#") == NULL) { + sscanf(line,"%s %lf %lf %lf", FamName, &dx, &dy, &dr); + Fnum = ElemIndex(FamName); + if (Fnum > 0) { + SetTol(Fnum, dx, dy, dr); + } else { + printf("LoadTolerances: undefined element %s\n", FamName); + exit_(1); + } + } + } while (fgets(line, 128, tolfile) != NULL); + + fclose(tolfile); +} + + +// Load tolerances from the file +void ScaleTolerances(const char *TolFileName, const double scl) +{ + char line[128], FamName[32]; + int Fnum; + double dx, dy, dr; + FILE *tolfile; + + tolfile = file_read(TolFileName); + + do + fgets(line, 128, tolfile); + while (strstr(line, "#") != NULL); + + do { + if (strstr(line, "#") == NULL) { + sscanf(line,"%s %lf %lf %lf", FamName, &dx, &dy, &dr); + Fnum = ElemIndex(FamName); + if (Fnum > 0) { + Scale_Tol(Fnum, scl*dx, scl*dy, scl*dr); + } else { + printf("ScaleTolerances: undefined element %s\n", FamName); + exit_(1); + } + } + } while (fgets(line, 128, tolfile) != NULL); + fclose(tolfile); +} + + +void SetKpar(int Fnum, int Knum, int Order, double k) +{ + + Cell[Elem_GetPos(Fnum, Knum)].Elem.M->PBpar[Order+HOMmax] = k; + Mpole_SetPB(Fnum, Knum, Order); +} + + +void SetL(int Fnum, int Knum, double L) +{ + + Cell[Elem_GetPos(Fnum, Knum)].Elem.PL = L; +} + + +void SetL(int Fnum, double L) +{ + int i; + + for (i = 1; i <= GetnKid(Fnum); i++) + Cell[Elem_GetPos(Fnum, i)].Elem.PL = L; +} + + +void SetdKpar(int Fnum, int Knum, int Order, double dk) +{ + + Cell[Elem_GetPos(Fnum, Knum)].Elem.M->PBpar[Order+HOMmax] += dk; + Mpole_SetPB(Fnum, Knum, Order); +} + + +void SetKLpar(int Fnum, int Knum, int Order, double kL) +{ + long int loc; + + loc = Elem_GetPos(Fnum, Knum); + if (Cell[loc].Elem.PL != 0e0) + Cell[loc].Elem.M->PBpar[Order+HOMmax] = kL/Cell[loc].Elem.PL; + else + Cell[loc].Elem.M->PBpar[Order+HOMmax] = kL; + Mpole_SetPB(Fnum, Knum, Order); +} + + +void SetdKLpar(int Fnum, int Knum, int Order, double dkL) +{ + long int loc; + + loc = Elem_GetPos(Fnum, Knum); + if (Cell[loc].Elem.PL != 0e0) + Cell[loc].Elem.M->PBpar[Order + HOMmax] += dkL/Cell[loc].Elem.PL; + else + Cell[loc].Elem.M->PBpar[Order + HOMmax] += dkL; + Mpole_SetPB(Fnum, Knum, Order); +} + + +void SetdKrpar(int Fnum, int Knum, int Order, double dkrel) +{ + long int loc; + + loc = Elem_GetPos(Fnum, Knum); + if (Order == Dip && Cell[loc].Elem.M->Pthick == thick) + Cell[loc].Elem.M->PBpar[Dip+HOMmax] += dkrel*Cell[loc].Elem.M->Pirho; + else + Cell[loc].Elem.M->PBpar[Order+HOMmax] + += dkrel*Cell[loc].Elem.M->PBpar[Order+HOMmax]; + Mpole_SetPB(Fnum, Knum, Order); +} + + +void Setbn(int Fnum, int order, double bn) +{ + int i; + + for (i = 1; i <= GetnKid(Fnum); i++) + SetKpar(Fnum, i, order, bn); +} + + +void SetbnL(int Fnum, int order, double bnL) +{ + int i; + + for (i = 1; i <= GetnKid(Fnum); i++) + SetKLpar(Fnum, i, order, bnL); +} + + +void Setdbn(int Fnum, int order, double dbn) +{ + int i; + + for (i = 1; i <= GetnKid(Fnum); i++) + SetdKpar(Fnum, i, order, dbn); +} + + +void SetdbnL(int Fnum, int order, double dbnL) +{ + int i; + + for (i = 1; i <= GetnKid(Fnum); i++) { + SetdKLpar(Fnum, i, order, dbnL); + } +} + + +void Setbnr(int Fnum, long order, double bnr) +{ + int i; + + for (i = 1; i <= GetnKid(Fnum); i++) + SetdKrpar(Fnum, i, order, bnr); +} + + +void SetbnL_sys(int Fnum, int Order, double bnL_sys) +{ + int Knum; + long int loc; + + for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) { + loc = Elem_GetPos(Fnum, Knum); + if (Cell[loc].Elem.PL != 0.0) + Cell[loc].Elem.M->PBsys[Order+HOMmax] = bnL_sys/Cell[loc].Elem.PL; + else + Cell[loc].Elem.M->PBsys[Order+HOMmax] = bnL_sys; + Mpole_SetPB(Fnum, Knum, Order); + } +} + + +void set_dbn_rel(const int type, const int n, const double dbn_rel) +{ + long int j; + double dbn; + + printf("\n"); + printf("Setting Db_%d/b_%d = %6.1e for:\n", n, type, dbn_rel); + printf("\n"); + for (j = 0; j <= globval.Cell_nLoc; j++) + if ((Cell[j].Elem.Pkind == Mpole) && (Cell[j].Elem.M->n_design == type)) { + printf("%s\n", Cell[j].Elem.PName); + dbn = dbn_rel*Cell[j].Elem.M->PBpar[type+HOMmax]; + Cell[j].Elem.M->PBrms[n+HOMmax] = dbn; + Cell[j].Elem.M->PBrnd[n+HOMmax] = normranf(); + Mpole_SetPB(Cell[j].Fnum, Cell[j].Knum, n); + } +} + + +double GetKpar(int Fnum, int Knum, int Order) +{ + return (Cell[Elem_GetPos(Fnum, Knum)].Elem.M->PBpar[Order+HOMmax]); +} + + +double GetL(int Fnum, int Knum) +{ + return (Cell[Elem_GetPos(Fnum, Knum)].Elem.PL); +} + + +double GetKLpar(int Fnum, int Knum, int Order) +{ + long int loc; + + loc = Elem_GetPos(Fnum, Knum); + if (Cell[loc].Elem.PL != 0e0) + return (Cell[loc].Elem.M->PBpar[Order+HOMmax]*Cell[loc].Elem.PL); + else + return (Cell[loc].Elem.M->PBpar[Order+HOMmax]); +} + + +void SetdKLrms(int Fnum, int Order, double dkLrms) +{ + long int Knum, loc; + + for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) { + loc = Elem_GetPos(Fnum, Knum); + if (Cell[loc].Elem.PL != 0e0) + Cell[loc].Elem.M->PBrms[Order+HOMmax] = dkLrms/Cell[loc].Elem.PL; + else + Cell[loc].Elem.M->PBrms[Order+HOMmax] = dkLrms; + Cell[loc].Elem.M->PBrnd[Order+HOMmax] = normranf(); + Mpole_SetPB(Fnum, Knum, Order); + } +} + + +void Setdkrrms(int Fnum, int Order, double dkrrms) +{ + long int Knum, loc; + + for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) { + loc = Elem_GetPos(Fnum, Knum); + if (Order == Dip && Cell[loc].Elem.M->Pthick == thick) + Cell[loc].Elem.M->PBrms[Dip+HOMmax] = dkrrms*Cell[loc].Elem.M->Pirho; + else + Cell[loc].Elem.M->PBrms[Order+HOMmax] + = dkrrms*Cell[loc].Elem.M->PBpar[Order+HOMmax]; + Cell[loc].Elem.M->PBrnd[Order+HOMmax] = normranf(); + Mpole_SetPB(Fnum, Knum, Order); + } +} + + +void SetKL(int Fnum, int Order) +{ + long int Knum; + + for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) + Mpole_SetPB(Fnum, Knum, Order); +} + + +void set_dx(const int type, const double sigma_x, const double sigma_y) +{ + long int j; + + printf("\n"); + printf("Setting sigma_x,y = (%6.1e, %6.1e) for b_%d:\n", + sigma_x, sigma_y, type); + printf("\n"); + for (j = 0; j <= globval.Cell_nLoc; j++) + if ((Cell[j].Elem.Pkind == Mpole) && (Cell[j].Elem.M->n_design == type)) { + printf("%s\n", Cell[j].Elem.PName); + Cell[j].Elem.M->PdSrms[X_] = sigma_x; + Cell[j].Elem.M->PdSrms[Y_] = sigma_y; + Cell[j].Elem.M->PdSrnd[X_] = normranf(); + Cell[j].Elem.M->PdSrnd[Y_] = normranf(); + Mpole_SetdS(Cell[j].Fnum, Cell[j].Knum); + } +} + + +void SetBpmdS(int Fnum, double dxrms, double dyrms) +{ + long int Knum, loc; + + for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) { + loc = Elem_GetPos(Fnum, Knum); + Cell[loc].dS[X_] = normranf()*dxrms; Cell[loc].dS[Y_] = normranf()*dyrms; + } +} + + +/****************************************/ +/* Routines for closed orbit correction */ +/****************************************/ + + +/****************************************************************************/ +void codstat(double *mean, double *sigma, double *xmax, long lastpos, bool all) +{ + long i, j, n; + Vector2 sum, sum2; + double TEMP; + + n = 0; + for (j = 0; j <= 1; j++) { + sum[j] = 0.0; sum2[j] = 0.0; xmax[j] = 0.0; + } + for (i = 0; i <= lastpos; i++) { + if (all || Cell[i].Fnum == globval.bpm) { + n++; + for (j = 1; j <= 2; j++) { + sum[j - 1] += Cell[i].BeamPos[j * 2 - 2]; + TEMP = Cell[i].BeamPos[j * 2 - 2]; + sum2[j - 1] += TEMP * TEMP; + xmax[j - 1] = max(xmax[j - 1], fabs(Cell[i].BeamPos[j * 2 - 2])); + } + } + } + for (j = 0; j <= 1; j++) { + if (n != 0) + mean[j] = sum[j] / n; + else + mean[j] = 0.0; + if (n != 0 && n != 1) { + TEMP = sum[j]; + sigma[j] = (n * sum2[j] - TEMP * TEMP) / (n * (n - 1.0)); + } else + sigma[j] = 0.0; + if (sigma[j] >= 0.0) + sigma[j] = sqrt(sigma[j]); + else + sigma[j] = 0.0; + } +} + +/****************************************************************************/ +/* void CodStatBpm(double *mean, double *sigma, double *xmax, long lastpos, + long bpmdis[mnp]) + + Purpose: + Get statistics for closed orbit + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + none + +****************************************************************************/ +void CodStatBpm(double *mean, double *sigma, double *xmax, long lastpos, + long bpmdis[mnp]) +{ + long i, j, m, n; + Vector2 sum, sum2; + double TEMP; + + m= n= 0; + for (j = 0; j <= 1; j++) { + sum[j] = 0.0; sum2[j] = 0.0; xmax[j] = 0.0; + } + + for (i = 0; i <= lastpos; i++) { + if (Cell[i].Fnum == globval.bpm) { + if (! bpmdis[m]) { + for (j = 1; j <= 2; j++) { + sum[j - 1] += Cell[i].BeamPos[j * 2 - 2]; + TEMP = Cell[i].BeamPos[j * 2 - 2]; + sum2[j - 1] += TEMP * TEMP; + xmax[j - 1] = max(xmax[j - 1], fabs(Cell[i].BeamPos[j * 2 - 2])); + } + n++; + } + m++; + } + } + for (j = 0; j <= 1; j++) { + if (n != 0) + mean[j] = sum[j] / n; + else + mean[j] = 0.0; + if (n != 0 && n != 1) { + TEMP = sum[j]; + sigma[j] = (n * sum2[j] - TEMP * TEMP) / (n * (n - 1.0)); + } else + sigma[j] = 0.0; + if (sigma[j] >= 0.0) + sigma[j] = sqrt(sigma[j]); + else + sigma[j] = 0.0; + } +} + + + +/****************************************************************************/ +/* double digitize(double x, double maxkick, double maxsamp) + + Purpose: + Map x onto the integer interval (-maxsamp ... maxsamp) where maxsamp + corresponds maxkick. + + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + none + +****************************************************************************/ +double digitize(double x, double maxkick, double maxsamp) +{ + if (maxkick>0.) + if (maxsamp>1.) + return Sgn(x)*maxkick/maxsamp *min(floor(fabs(x)/maxkick*maxsamp),maxsamp-1.); + else { + return Sgn(x)*min(fabs(x),maxkick); + } + else + return x; +} + + +/****************************************************************************/ +/* double digitize2(long plane, long inum, double x, double maxkick, double maxsamp) + + Purpose: + + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + none + +****************************************************************************/ +svdarray xmemo[2]; + +double digitize2(long plane, long inum, double x, double maxkick, double maxsamp) +{ + double xint; + + if (maxkick>0.) + if (maxsamp>1.) + { + xint=min(floor(fabs(x)/maxkick*maxsamp),maxsamp-1.); + + if(fabs(xint-xmemo[inum][plane]) >=1.) + { + xmemo[inum][plane]=xint; + } + else + { + xmemo[inum][plane]+=0.1; + xint=min(xmemo[inum][plane],maxsamp-1.); + } + return Sgn(x)*maxkick/maxsamp*xint; + } + else + { + return Sgn(x)*min(fabs(x),maxkick); + } + else + return x; +} + + +// MATH ROUTINE a mettre dans mathlib.c + +/****************************************************************************/ +/* void GetMean(n, x) + + Purpose: + Get out the mean value of vector x + + Input: + n vector size + x vector to get out the mean value + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + to be moved in mathlib + +****************************************************************************/ +void GetMean(long n, double *x) +{ + long i; + double mean = 0e0; + + if ( n < 1 ) + { + fprintf(stdout,"GetMean: error wrong vector size n=%ld\n",n); + exit_(1); + } + for (i = 0; i < n; i++) + mean += x[i]; + mean /= n; + for (i = 0; i < n; i++) + x[i] = x[i] - mean; +} + +/****************************************************************************/ +/* double Fract(double x) + + Purpose: + Gets fractional part of x + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + none + +****************************************************************************/ +double Fract(double x) +{ + return (x - (long int)x); +} + +/****************************************************************************/ +/* double Sgn (double x) + + Purpose: + Gets sign of x + + Input: + none + + Output: + 0 if zero + 1 if positive + -1 if negative + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + none + +****************************************************************************/ +double Sgn (double x) +{ + return (x == 0.0 ? 0.0 : (x > 0.0 ? 1.0 : -1.0)); +} + + +void ChamberOff(void) +{ + int i; + + for (i = 0; i <= globval.Cell_nLoc; i++) { + Cell[i].maxampl[X_][0] = -max_ampl; + Cell[i].maxampl[X_][1] = max_ampl; + Cell[i].maxampl[Y_][0] = -max_ampl; + Cell[i].maxampl[Y_][1] = max_ampl; + } + status.chambre = false; +} + + + void PrintCh(void) +{ + long i = 0; + struct tm *newtime; + FILE *f; + + const char *fic = "chambre.out"; + + newtime = GetTime(); + + f = file_write(fic); + //fprintf(f, "# TRACY III synchrotron soleil -- %s -- %s \n", fic, asctime2(newtime)); + //fprintf(f, "# name s -xch +xch zch\n"); + //fprintf(f, "# [mm] [mm] [mm]\n"); + //fprintf(f, "#\n"); + fprintf(f, "# i name s -xch(mm) +xch(mm) zch (mm)\n#\n"); + + for (i = 0; i <= globval.Cell_nLoc; i++) + // fprintf(f, "%4ld %15s %6.2f %7.3f %7.3f %7.3f\n", + // i, Cell[i].Elem.PName, Cell[i].S, +// Cell[i].maxampl[X_][0]*1E3, Cell[i].maxampl[X_][1]*1E3, +// Cell[i].maxampl[Y_][1]*1E3); + + fprintf(f, "%4ld ", i); + fprintf(f, "%6.2f %7.3f %7.3f %7.3f\n", Cell[i].S, + Cell[i].maxampl[X_][0] * 1E3, + Cell[i].maxampl[X_][1] * 1E3, + Cell[i].maxampl[Y_][1] * 1E3); + + + fclose(f); +} + + + +// /** function from soleilcommon.c **/ +// Version of nsrl-ii +// Move Read_Lattice() to soleilcommon.cc +//void Read_Lattice(char *fic) +//{ +// bool status; +// char fic_maille[S_SIZE+4] = "", fic_erreur[S_SIZE+4] = ""; +// int i; +// double dP = 0.0; +// Vector2 beta, alpha, eta, etap; +// Vector codvect; + +// const double RFacceptance = 0.060001; // soleil energy acceptance + +// strcpy(fic_maille, fic); strcpy(fic_erreur, fic); + +// /* generation automatique du nom du fichier maille et erreur */ +// strcat(fic_maille, ".lat"); strcat(fic_erreur, ".lax"); + +// /* Initialisation de Tracy */ + +// t2init(); + + /* open the lattice Input file */ + +// if ((fi = fopen(fic_maille, "r")) == NULL) { +// fprintf(stdout, "ReadLattice: Error while opening file %s \n", fic_maille); +// fprintf(stdout, "The lattice file has to end by .lat \n"); +// exit_(1); +// } + +// /* opens the lattice Output file */ +// // if ((fo = fopen(fic_erreur, "w")) == NULL) { +// // fprintf(stdout, "ReadLattice: Error while opening file %s \n", fic_erreur); +// exit_(1); +// } + + /* Reads lattice and set principle parameters + * Energy CODeps and energy offset + * print statistics + */ +// status = Lattice_Read(&fi, &fo); + +// if (status == false) { +// cout << "Lattice_Read function has returned false" << endl; +// cout << "See file " << fic_erreur << endl; +// exit_(1); +// } +// cout << "Lattice file: " << fic_maille << endl; + + /* initializes cell structure: construction of the RING */ + /* Creator of all the matrices for each element */ +// Cell_Init(); + +// if (globval.RingType == 1) { // for a ring +// /* define x/y physical aperture */ +// ChamberOff(); + +// /* Defines global variables for Tracy code */ +// globval.H_exact = false; // Small Ring Hamiltonian +// globval.quad_fringe = false; // quadrupole fringe fields on/off +// globval.EPU = false; // Elliptically Polarizing Undulator +// globval.Cavity_on = false; /* Cavity on/off */ +// globval.radiation = false; /* radiation on/off */ +// globval.IBS = false; /* diffusion on/off */ +// globval.emittance = false; /* emittance on/off */ +// globval.pathlength = false; /* Path lengthening computation */ +// globval.CODimax = 40; /* maximum number of iterations for COD +// algo */ +// globval.delta_RF = RFacceptance;/* energy acceptance for SOLEIL */ + +// Cell_SetdP(dP); /* added for correcting BUG if non convergence: +// compute on momentum linear matrices */ +// } else { + // for transfer lines + /* Initial settings : */ +// beta[X_] = 8.1; +// alpha[X_] = 0.0; +// beta[Y_] = 8.1; +// alpha[Y_] = 0.0; +// eta[X_] = 0.0; +// etap[X_] = 0.0; +// eta[Y_] = 0.0; +// etap[Y_] = 0.0; + +// for (i = 0; i < ss_dim; i++) { +// codvect[i] = 0.0; +// globval.CODvect[i] = codvect[i]; +// } +// dP = codvect[delta_]; + +// /* Defines global variables for Tracy code */ +// globval.Cavity_on = false; /* Cavity on/off */ +// globval.radiation = false; /* radiation on/off */ +// globval.emittance = false; /* emittance on/off */ +// globval.pathlength = false; /* Path lengthening computation */ + // globval.CODimax = 10; /* maximum number of iterations for COD +// algo */ +// globval.delta_RF = RFacceptance; /* 6% + epsilon energy acceptance +// for SOLEIL */ +// globval.dPparticle = dP; + +// ChamberOff(); + +// TransTwiss(alpha, beta, eta, etap, codvect); +// } +//} + + + +/****************************************************************************/ +/* void GetChromTrac(long Nb, long Nbtour, double emax, + double *xix, double *xiz) + + Purpose: + Computes chromaticities by tracking + + Input: + Nb point number + Nbtour turn number + emax energy step + + Output: + xix horizontal chromaticity + xiz vertical chromaticity + + Return: + none + + Global variables: + trace + + Specific functions: + Trac_Simple, Get_NAFF + + Comments: + 27/04/03 chromaticities are now output arguments + +****************************************************************************/ +#define nterm 2 +void GetChromTrac(long Nb, long Nbtour, double emax, double *xix, double *xiz) +{ + bool status = true; + int nb_freq[2] = { 0, 0 }; /* frequency number to look for */ + int i = 0; + double Tab[6][NTURN], fx[nterm], fz[nterm], nux1, nux2, nuz1, nuz2; + + double x = 1e-6, xp = 0.0, z = 1e-6, zp = 0.0; + double x0 = 1e-6, xp0 = 0.0, z0 = 1e-6, zp0 = 0.0; + + /* initializations */ + for (i = 0; i < nterm; i++) { + fx[i] = 0.0; fz[i] = 0.0; + } + /* end init */ + + /* Tracking for delta = emax and computing tunes */ + x = x0; xp = xp0; z = z0; zp = zp0; + + Trac_Simple(x, xp, z, zp, emax, 0.0, Nbtour, Tab, &status); + Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq); + + nux1 = (fabs (fx[0]) > 1e-8 ? fx[0] : fx[1]); nuz1 = fz[0]; + + if (trace) + fprintf(stdout, + "\n Entering routine for chroma using tracking\n"); + if (trace) + fprintf(stdout, "emax= % 10.6e nux1=% 10.6e nuz1= % 10.6e\n", + emax, nux1, nuz1); + + /* Tracking for delta = -emax and computing tunes */ + x = x0; xp = xp0; z = z0; zp = zp0; + + Trac_Simple(x, xp, z, zp, -emax, 0.0, Nbtour, Tab, &status); + Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq); + + if (trace) + fprintf(stdout, "nturn=%6ld x=% 10.5g xp=% 10.5g z=% 10.5g zp=% 10.5g" + " delta=% 10.5g ctau=% 10.5g \n", + Nbtour, + Tab[0][Nbtour-1], Tab[1][Nbtour-1], + Tab[2][Nbtour-1], Tab[3][Nbtour-1], + Tab[4][Nbtour-1], Tab[5][Nbtour-1]); + + nux2 = (fabs(fx[0]) > 1e-8 ? fx[0] : fx[1]); nuz2 = fz[0]; + + if (trace) + fprintf(stdout, "emax= % 10.6e nux2= % 10.6e nuz2= % 10.6e\n", + -emax, nux2, nuz2); + + /* Computing chromaticities */ + *xix = (nux2-nux1)*0.5/emax; *xiz = (nuz2-nuz1)*0.5/emax; + + if (trace) + fprintf (stdout, " Exiting routine for chroma using tracking\n\n"); +} +#undef nterm + +/****************************************************************************/ +/* void GetTuneTrac(long Nbtour, double emax, double *nux, double *nuz) + + Purpose: + Computes chromaticities by tracking + + Input: + Nb point number + Nbtour turn number + emax energy step + + Output: + none + + Return: + none + + Global variables: + trace + + Specific functions: + Trac_Simple, Get_NAFF + + Comments: + none + +****************************************************************************/ +#define nterm 2 +void GetTuneTrac(long Nbtour, double emax, double *nux, double *nuz) +{ + double Tab[6][NTURN], fx[nterm], fz[nterm]; + int nb_freq[2]; + bool status; + + double x = 1e-6, xp = 0.0, z = 1e-6, zp = 0.0; + + Trac_Simple(x, xp, z, zp, emax, 0.0, Nbtour, Tab, &status); + Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq); + + *nux = (fabs (fx[0]) > 1e-8 ? fx[0] : fx[1]); + *nuz = fz[0]; +} +#undef nterm + +/****************************************************************************/ +/* void TransTwiss(double *alpha, double *beta, double *eta, double *etap, double *codvect) + + Purpose: high level application + Calculate Twiss functions for a transport line + + Input: + alpha alpha fonctions at the line entrance + beta beta fonctions at the line entrance + eta disperion fonctions at the line entrance + etap dipersion derivatives fonctions at the line entrance + codvect closed orbit fonctions at the line entrance + + Output: + none + + Return: + none + + Global variables: + + + Specific functions: + TransTrace + + Comments: + redundant with ttwiss + +****************************************************************************/ +void TransTwiss(Vector2 &alpha, Vector2 &beta, Vector2 &eta, Vector2 &etap, + Vector &codvect) +{ + TransTrace(0, globval.Cell_nLoc, alpha, beta, eta, etap, codvect); +} + + +/****************************************************************************/ +/* void ttwiss(double *alpha, double *beta, double *eta, double *etap, double dP) + + Purpose: + Calculate Twiss functions for transport line + + Input: + none + + Output: + none + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + redundant with TransTwiss + +****************************************************************************/ +void ttwiss(const Vector2 &alpha, const Vector2 &beta, + const Vector2 &eta, const Vector2 &etap, const double dP) +{ + TraceABN(0, globval.Cell_nLoc, alpha, beta, eta, etap, dP); +} + +/****************************************************************************/ +/* void findcodS(double dP) + + Purpose: + Search for the closed orbit using a numerical method + Algo: Newton_Raphson method + Quadratic convergence + May need a guess starting point + Simple precision algorithm + + Input: + dP energy offset + + Output: + none + + Return: + none + + Global variables: + none + + specific functions: + Newton_Raphson + + Comments: + Method introduced because of bad convergence of da for ID using RADIA maps + +****************************************************************************/ +void findcodS(double dP) +{ + double *vcod; + Vector x0; + const int ntrial = 40; // maximum number of trials for closed orbit + const double tolx = 1e-8; // numerical precision + int k; + int dim; // 4D or 6D tracking + long lastpos; + + vcod = dvector(1, 6); + + // starting point + for (k = 1; k <= 6; k++) + vcod[k] = 0.0; + + vcod[5] = dP; // energy offset + + if (globval.Cavity_on){ + dim = 6; /* 6D tracking*/ + fprintf(stdout,"Error looking for cod in 6D\n"); + exit_(1); + } + else{ + dim = 4; /* 4D tracking */ + vcod[1] = Cell[0].Eta[0]*dP; vcod[2] = Cell[0].Etap[0]*dP; + vcod[3] = Cell[0].Eta[1]*dP; vcod[4] = Cell[0].Etap[1]*dP; + } + + Newton_RaphsonS(ntrial, vcod, dim, tolx); + + if (status.codflag == false) + fprintf(stdout, "Error No COD found\n"); + if (trace) { + for (k = 1; k <= 6; k++) + x0[k-1] = vcod[k]; + fprintf(stdout, "Before cod % .5e % .5e % .5e % .5e % .5e % .5e \n", + x0[0], x0[1], x0[2], x0[3], x0[4], x0[5]); + Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); + fprintf(stdout, "After cod % .5e % .5e % .5e % .5e % .5e % .5e \n", + x0[0], x0[1], x0[2], x0[3], x0[4], x0[5]); + Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); + } + free_dvector(vcod,1,6); +} + +/****************************************************************************/ +/* void findcod(double dP) + + Purpose: + Search for the closed orbit using a numerical method + Algo: Newton_Raphson method + Quadratic convergence + May need a guess starting point + Simple precision algorithm + 4D + Starting point: linear closed orbit + + 6D + Starting point: zero + if radiation on : x[5] is the synchroneous phase (equilibrium RF phase) + off: x[5] is zero + + Input: + dP energy offset + + Output: + none + + Return: + none + + Global variables: + none + + specific functions: + Newton_Raphson + + Comments: + Method introduced because of bad convergence of da for ID + using RADIA maps + +****************************************************************************/ +void findcod(double dP) +{ + Vector vcod; + const int ntrial = 40; // maximum number of trials for closed orbit + const double tolx = 1e-10; // numerical precision + int k, dim = 0; + long lastpos; + + // initializations + for (k = 0; k <= 5; k++) + vcod[k] = 0.0; + + if (globval.Cavity_on){ + fprintf(stdout,"warning looking for cod in 6D\n"); + dim = 6; + } else{ // starting point linear closed orbit + dim = 4; + vcod[0] = Cell[0].Eta[0]*dP; vcod[1] = Cell[0].Etap[0]*dP; + vcod[2] = Cell[0].Eta[1]*dP; vcod[3] = Cell[0].Etap[1]*dP; + vcod[4] = dP; // energy offset + } + + Newton_Raphson(dim, vcod, ntrial, tolx); + + if (status.codflag == false) + fprintf(stdout, "Error No COD found\n"); + + CopyVec(6, vcod, globval.CODvect); // save closed orbit at the ring entrance + + if (trace) + { + fprintf(stdout, + "Before cod2 % .5e % .5e % .5e % .5e % .5e % .5e \n", + vcod[0], vcod[1], vcod[2], vcod[3], vcod[4], vcod[5]); + Cell_Pass(0, globval.Cell_nLoc, vcod, lastpos); + fprintf(stdout, + "After cod2 % .5e % .5e % .5e % .5e % .5e % .5e \n", + vcod[0], vcod[1], vcod[2], vcod[3], vcod[4], vcod[5]); + } +} +/****************************************************************************/ +/* void computeFandJS(double *x, int n, double **fjac, double *fvect) + + Purpose: + Simple precision algo + Tracks x over one turn. And computes the Jacobian matrix of the + transformation by numerical differentiation. + using forward difference formula : faster but less accurate + using symmetric difference formula + + Input: + x vector for evaluation + n dimension 4 or 6 + + Output: + fvect transport of x over one turn + fjac Associated jacobian matrix + + Return: + none + + Global variables: + none + + specific functions: + none + + Comments: + none + +****************************************************************************/ + +void computeFandJS(double *x, int n, double **fjac, double *fvect) +{ + int i, k; + long lastpos = 0L; + Vector x0, fx, fx1, fx2; + + const double deps = 1e-8; //stepsize for numerical differentiation + + for (i = 1; i <= 6; i++) + x0[i - 1] = x[i]; + + Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); + + for (i = 1; i <= n; i++) + { + fvect[i] = x0[i - 1]; + fx[i - 1] = x0[i - 1]; + } + + // compute Jacobian matrix by numerical differentiation + for (k = 0; k < n; k++) + { + for (i = 1; i <= 6; i++) + x0[i - 1] = x[i]; + x0[k] += deps; // differential step in coordinate k + + Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); // tracking along the ring + for (i = 1; i <= 6; i++) + fx1[i - 1] = x0[i - 1]; + + for (i = 1; i <= 6; i++) + x0[i - 1] = x[i]; + x0[5] = 0.0; + x0[k] -= deps; // differential step in coordinate k + + Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); // tracking along the ring + for (i = 1; i <= 6; i++) + fx2[i - 1] = x0[i - 1]; + + for (i = 1; i <= n; i++) // symmetric difference formula + fjac[i][k + 1] = 0.5 * (fx1[i - 1] - fx2[i - 1]) / deps; + //~ for (i = 1; i <= n; i++) // forward difference formula + //~ fjac[i][k + 1] = (float) ((x0[i - 1] - fx[i - 1]) / deps); + } +} + +/****************************************************************************/ +/* void computeFand(int n, float *x, float **fjac, float *fvect) + + Purpose: + Tracks x over one turn. And computes the Jacobian matrix of the + transformation by numerical differentiation. + using symmetric difference formula + double precision algorithm + + Input: + x vector for evaluation + + Output: + fvect transport of x over one turn + fjac Associated jacobian matrix + + Return: + none + + Global variables: + none + + specific functions: + none + + Comments: + none + +****************************************************************************/ +void computeFandJ(int n, Vector &x, Matrix &fjac, Vector &fvect) +{ + int i, k; + long lastpos = 0; + Vector x0, fx1, fx2; + + const double deps = 1e-8; //stepsize for numerical differentiation + + CopyVec(6, x, x0); + + Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); + CopyVec(n, x0, fvect); + + // compute Jacobian matrix by numerical differentiation + for (k = 0; k < n; k++) { + CopyVec(6L, x, x0); + x0[k] += deps; // differential step in coordinate k + + Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); // tracking along the ring + CopyVec(6L, x0, fx1); + + CopyVec(6L, x, x0); + x0[k] -= deps; // differential step in coordinate k + + Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); // tracking along the ring + CopyVec(6L, x0, fx2); + + for (i = 0; i < n; i++) // symmetric difference formula + fjac[i][k] = 0.5 * (fx1[i] - fx2[i]) / deps; + } +} + +/****************************************************************************/ +/* void Newton_RaphsonS(int ntrial,double x[],int n,double tolx, double tolf) + + Purpose: + Newton_Rapson algorithm from Numerical Recipes + single precision algorithm + Robustess: quadratic convergence + Hint: for n-dimensional problem, the algo can be stuck on local minimum + In this case, it should be enough to provide a resonable starting + point. + + Method: + look for closed orbit solution of f(x) = x + This problems is equivalent to finding the zero of g(x)= f(x) - x + g(x+h) ~= f(x) - x + (Jacobian(f) -Id) h + O(h*h) + Then at first order we solve h: + h = - inverse(Jacobian(f) -Id) * (f(x)-x) + the new guess is then xnew = x + h + By iteration, this converges quadratically. + + The algo is stopped whenever |x -xnew| < tolx + + f(x) is computes by tracking over one turn + Jacobian(f) is computed numerically by numerical differentiation + These two operations are provided by the function computeFandJ + + Input: + ntrial number of iterations for closed zero search + n number of dimension 4 or 6 + x intial guess for the closed orbit + tolx tolerance over the solution x + tolf tolerance over the evalution f(x) + + Output: + x closed orbit + + Return: + none + + Global variables: + status + + specific functions: + computeFandJS + ludcmp,lubksb + + Comments: + none + +****************************************************************************/ + +void Newton_RaphsonS(int ntrial, double x[], int n, double tolx) +{ + int k, i, *indx; + double errx, d, *bet, *fvect, **alpha; + + errx = 0.0; + // NR arrays start from 1 and not 0 !!! + indx = ivector(1, n); + bet = dvector(1, n); + fvect = dvector(1, n); + alpha = dmatrix(1, n, 1, n); + + for (k = 1; k <= ntrial; k++) { // loop over number of iterations + // supply function values at x in fvect and Jacobian matrix in fjac + computeFandJS(x, n, alpha, fvect); + + // Jacobian -Id + for (i = 1; i <= n; i++) + alpha[i][i] -= 1.0; + for (i = 1; i <= n; i++) + bet[i] = x[i] - fvect[i]; // right side of linear equation + // solve linear equations using LU decomposition using NR routines + dludcmp(alpha, n, indx, &d); + dlubksb(alpha, n, indx, bet); + errx = 0.0; // check root convergence + for (i = 1; i <= n; i++) { // update solution + errx += fabs(bet[i]); + x[i] += bet[i]; + } + + if (trace) + fprintf(stdout, + "%02d: cod % .5e % .5e % .5e % .5e % .5e % .5e errx =% .5e\n", + k, x[1], x[2], x[3], x[4], x[5], x[6], errx); + if (errx <= tolx) { + status.codflag = true; + break; + } + } + // check whever closed orbit found out + if ((k >= ntrial) && (errx >= tolx * 100)) status.codflag = false; + + free_dmatrix(alpha,1,n,1,n); free_dvector(bet,1,n); free_dvector(fvect,1,n); + free_ivector(indx,1,n); +} + + +/****************************************************************************/ +/* int Newton_Raphson(int n, double x[], int ntrial, double tolx) + + Purpose: + Newton_Rapson algorithm from Numerical Recipes + double precision algorithm + Robustess: quadratic convergence + Hint: for n-dimensional problem, the algo can be stuck on local minimum + In this case, it should be enough to provide a resonable starting + point. + + Method: + look for closed orbit solution of f(x) = x + This problems is equivalent to finding the zero of g(x)= f(x) - x + g(x+h) ~= f(x) - x + (Jacobian(f) -Id) h + O(h*h) + Then at first order we solve h: + h = - inverse(Jacobian(f) -Id) * (f(x)-x) + the new guess is then xnew = x + h + By iteration, this converges quadratically. + + The algo is stopped whenever |x -xnew| < tolx + + f(x) is computes by tracking over one turn + Jacobian(f) is computed numerically by numerical differentiation + These two operations are provided by the function computeFandJ + + Input: + ntrial number of iterations for closed zero search + x intial guess for the closed orbit + tolx tolerance over the solution x + tolf tolerance over the evalution f(x) + + Output: + x closed orbit + + Return: + none + + Global variables: + status + + specific functions: + computeFandJ + InvMat, LinTrans + + Comments: + none + +****************************************************************************/ +int Newton_Raphson (int n, Vector &x, int ntrial, double tolx) +{ + int k, i; + double errx; + Vector bet, fvect; + Matrix alpha; + + errx = 0.0; + + for (k = 1; k <= ntrial; k++) { // loop over number of iterations + // supply function values at x in fvect and Jacobian matrix in fjac + computeFandJ(n, x, alpha, fvect); + + // Jacobian - Id + for (i = 0; i < n; i++) + alpha[i][i] -= 1.0; + for (i = 0; i < n; i++) + bet[i] = x[i] - fvect[i]; // right side of linear equation + // inverse matrix using gauss jordan method from Tracy (from NR) + if (!InvMat((long) n,alpha)) + fprintf(stdout,"Matrix non inversible ...\n"); + LinTrans((long) n, alpha, bet); // bet = alpha*bet + errx = 0.0; // check root convergence + for (i = 0; i < n; i++) + { // update solution + errx += fabs(bet[i]); + x[i] += bet[i]; + } + + if (trace) + fprintf(stdout, + "%02d: cod2 % .5e % .5e % .5e % .5e % .5e % .5e errx =% .5e\n", + k, x[0], x[1], x[2], x[3], x[4], x[5], errx); + if (errx <= tolx) + { + status.codflag = true; + return 1; + } + } + // check whever closed orbit found out + if ((k >= ntrial) && (errx >= tolx)) + { + status.codflag = false; + return 1; + } + return 0; +} + + + +void rm_mean(long int n, double x[]) +{ + long int i; + double mean; + + mean = 0.0; + for (i = 0; i < n; i++) + mean += x[i]; + mean /= n; + for (i = 0; i < n; i++) + x[i] -= mean; +} diff --git a/tracy/tracy/src/soleillib.cc b/tracy/tracy/src/soleillib.cc index 0b6ab72..1d57878 100644 --- a/tracy/tracy/src/soleillib.cc +++ b/tracy/tracy/src/soleillib.cc @@ -109,13 +109,16 @@ void InducedAmplitude(long spos) /* Ouverture fichier moustache */ if ((outf = fopen(nomfic, "w")) == NULL) { - fprintf(stdout, "Erreur � l'ouverture de %s\n",nomfic); + fprintf(stdout, "Error when open filename %s\n",nomfic); exit_(1); } + fprintf(outf, "# Induced amplitude transported at lattice entrance\n"); fprintf(outf, "# dp xind zind " - " Betax(0) Betaz(0) Betax betaz" - " Hx Hz etax etaxp\n#\n"); + " Betax(entrance) Betaz(entrance) Betax betaz" + " Hx(delta)/delta^2 Hz(delta)/delta^2 " + " Hx(delta) Hz(delta) etax(delta) etaxp(delta)\n#\n"); + lastpos = 1; @@ -147,14 +150,16 @@ void InducedAmplitude(long spos) amp[1] = codvector[spos][1]; fprintf(outf, "%+10.5e %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e " - "%+10.5e %+10.5e %+10.5e %+10.5e \n", - dP, codvector[spos][0], codvector[spos][1], Celldebut.Beta[0], Celldebut.Beta[1], Cell.Beta[0], Cell.Beta[1], H[0], H[1], Cell.Eta[0], Cell.Etap[0]); - } + "%+10.5e %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e \n", + dP, codvector[spos][0], codvector[spos][1], + Celldebut.Beta[0], Celldebut.Beta[1], Cell.Beta[0], Cell.Beta[1], + H[0], H[1], H[0]*dP20, H[1]*dP20, Cell.Eta[0], Cell.Etap[0]); + } fclose(outf); } /****************************************************************************/ -/* void Hfonction(long pos, double dP,Vector2 H) +/* void Hfonction(long pos, double dP) Purpose: Compute the Hfunction at position pos for the energy offset dP @@ -170,7 +175,7 @@ void InducedAmplitude(long spos) Indeed in general : xco = eta(dp)*dp + x0(defaults) Input: - none + pos: element index in the lattice. Output: none @@ -190,13 +195,15 @@ void InducedAmplitude(long spos) ****************************************************************************/ -void Hfonction(long pos, double dP,Vector2 H) +//void Hfonction(long pos, double dP,Vector2 H) +void Hfonction(long pos, double dP) { CellType Cell; long i; +Vector2 H; - Ring_GetTwiss(pos, dP); /* Compute and get Twiss parameters */ - getelem(pos, &Cell); /* Position sur l'element pos */ + Ring_GetTwiss(true, dP); /* Compute and get Twiss parameters */ + getelem(pos, &Cell); /* Position of the element pos */ i = 0; /* Horizontal */ H[i] = (1+Cell.Alpha[i]*Cell.Alpha[i])/Cell.Beta[i]*Cell.Eta[i]*Cell.Eta[i]+ @@ -209,7 +216,7 @@ void Hfonction(long pos, double dP,Vector2 H) } /****************************************************************************/ -/* void Hcofonction(long pos, double dP,Vector2 H) +/* void Hcofonction(long pos, double dP) Purpose: Compute the true Hfunction defined by the chromatic closed orbit @@ -239,18 +246,19 @@ void Hfonction(long pos, double dP,Vector2 H) Bug: Cell.BeamPos does not give closed orbit !!! ****************************************************************************/ -void Hcofonction(long pos, double dP,Vector2 H) +//void Hcofonction(long pos, double dP,Vector2 H) +void Hcofonction(long pos, double dP) { CellType Cell; long i; long lastpos = 1; - + Vector2 H; //~ getcod(dP, lastpos); /* determine closed orbit */ findcod(dP); if (lastpos != globval.Cell_nLoc) printf("Ring unstable for dp=%+e @ pos=%ld\n", dP, lastpos); - Ring_GetTwiss(pos, dP); /* Compute and get Twiss parameters */ - getelem(pos, &Cell); /* Position sur l'element pos */ + Ring_GetTwiss(true, dP); /* Compute and get Twiss parameters */ + getelem(pos, &Cell); /* Position of the element pos */ i = 0; /* Horizontal */ H[i] = (1+Cell.Alpha[i]*Cell.Alpha[i])/Cell.Beta[i]*Cell.BeamPos[i]*Cell.BeamPos[i]+ @@ -260,6 +268,7 @@ void Hcofonction(long pos, double dP,Vector2 H) H[i] = (1+Cell.Alpha[i]*Cell.Alpha[i])/Cell.Beta[i]*Cell.BeamPos[i+1]*Cell.BeamPos[i+1]+ 2*Cell.Alpha[i]*Cell.BeamPos[i+1]*Cell.BeamPos[i+2]+ Cell.Beta[i]*Cell.BeamPos[i+2]*Cell.BeamPos[i+2]; + fprintf(stdout,"H[0]=%10.6f,H[1]=%10.6f\n",H[0],H[1]); } @@ -386,37 +395,37 @@ void DefineCh(void) for (i = 0; i <= globval.Cell_nLoc; i++) { if (Cell[i].Elem.Pkind == marker){ if (strncmp(Cell[i].Elem.PName,"ssep",4) == 0){ - if (trace) fprintf(stdout,"trouve %s Element numero %ld \n", + if (trace) fprintf(stdout,"found %s Element number %ld \n", Cell[i].Elem.PName,i); isep1 = i; } if (strncmp(Cell[i].Elem.PName,"esep",4) == 0) { - if (trace) fprintf(stdout,"trouve %s Element numero %ld \n", + if (trace) fprintf(stdout,"found %s Element number %ld \n", Cell[i].Elem.PName,i-1); isep2 = i-1; } if (strncmp(Cell[i].Elem.PName,"ehu600",6) == 0) { - if (trace) fprintf(stdout,"trouve %s Element numero %ld \n", + if (trace) fprintf(stdout,"found %s Element number %ld \n", Cell[i].Elem.PName,i-1); hu600 = i-1; } if (strncmp(Cell[i].Elem.PName,"ssdm",4) == 0) { - if (trace) fprintf(stdout,"trouve %s Element numero %ld \n", + if (trace) fprintf(stdout,"found %s Element number %ld \n", Cell[i].Elem.PName,i); isdm1 = i; } if (strncmp(Cell[i].Elem.PName,"esdm",4) == 0) { - if (trace) fprintf(stdout,"trouve %s Element numero %ld \n", + if (trace) fprintf(stdout,"found %s Element number %ld \n", Cell[i].Elem.PName,i); isdm2 = i; } if (strncmp(Cell[i].Elem.PName,"ssdac",5) == 0) { - if (trace) fprintf(stdout,"trouve %s Element numero %ld \n", + if (trace) fprintf(stdout,"found %s Element number %ld \n", Cell[i].Elem.PName,i); isdac1 = i; } if (strncmp(Cell[i].Elem.PName,"esdac",5) == 0) { - if (trace) fprintf(stdout,"trouve %s Element numero %ld \n", + if (trace) fprintf(stdout,"found %s Element number %ld \n", Cell[i].Elem.PName,i-1); isdac2 = i-1; } @@ -433,23 +442,27 @@ void DefineCh(void) Cell[i].maxampl[X_][0] = -35.e-3; Cell[i].maxampl[X_][1] = 35.e-3; Cell[i].maxampl[Y_][1] = 12.5e-3; - } else if ((i >= isdm1) && (i <= isdm2)) { + } + else if ((i >= isdm1) && (i <= isdm2)) { /* SD13 */ Cell[i].maxampl[X_][0] = -21e-3; Cell[i].maxampl[X_][1] = 21e-3; Cell[i].maxampl[Y_][1] = 5e-3; - } else if ((i >= isep1) && (i <= isep2)) { + } + else if ((i >= isep1) && (i <= isep2)) { /* septum */ Cell[i].maxampl[X_][0] = -25e-3; Cell[i].maxampl[X_][1] = 25e-3; Cell[i].maxampl[Y_][1] = 12.5e-3; - } else if ((i >= isdac1) && (i <= isdac2)) { + } + else if ((i >= isdac1) && (i <= isdac2)) { /* minigap */ Cell[i].maxampl[X_][0] = -35e-3; Cell[i].maxampl[X_][1] = 35e-3; Cell[i].maxampl[Y_][1] = 2.5e-3; } - if (i <= hu600) { + /* if ((i>isep1) && (i<=hu600)) /* debut maille en section courte PB 25-09-07 */ + if (i <= hu600) { /* HU640 */ Cell[i].maxampl[Y_][1] = 7.0e-3; } @@ -482,11 +495,13 @@ void DefineCh(void) none ****************************************************************************/ +// Definition in soleilcommon.cc + void ChamberOn(void) { DefineCh(); status.chambre = true; -} +} /****************************************************************************/ /* void Trac_Tab(double x, double px, double y, double py, double dp, @@ -648,7 +663,7 @@ void NuDx(long Nbx, long Nbz, long Nbtour, double xmax, double zmax, exit_(1); } - fprintf(outf,"# Tracy-2 v. 2.8 -- %s -- %s \n", ficx, asctime2(newtime)); + fprintf(outf,"# Tracy III -- %s -- %s \n", ficx, asctime2(newtime)); fprintf(outf,"# nu = f(x) \n"); fprintf(outf,"# x[m] z[m] fx fz \n"); @@ -693,7 +708,7 @@ void NuDx(long Nbx, long Nbz, long Nbtour, double xmax, double zmax, exit_(1); } - fprintf(outf,"# tracy-3.5 -- %s -- %s \n", ficz, asctime2(newtime)); + fprintf(outf,"# tracy III -- %s -- %s \n", ficz, asctime2(newtime)); fprintf(outf,"# nu = f(z) \n"); fprintf(outf,"# x[mm] z[mm] fx fz \n"); @@ -741,13 +756,17 @@ double get_D(const double df_x, const double df_y) /****************************************************************************/ /* void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax, - double energy, bool diffusion) + double energy, bool diffusion, bool matlab) Purpose: Compute a frequency map of Nbx x Nbz points For each set of initial conditions the particle is tracked over Nbtour for an energy offset dp + Frequency map is based on fixed beam energy, trace x versus z, + or, tracking transverse dynamic aperture for fixed momentum + (usually, on-momentum) particle. + The stepsize follows a square root law Results in fmap.out @@ -759,7 +778,8 @@ double get_D(const double df_x, const double df_y) zmax vertical maximum amplitude Nbtour number of turn for tracking energy particle energy offset - + matlab set file print format for matlab post-process; specific for nsrl-ii + Output: status true if stable false otherwise @@ -779,8 +799,10 @@ double get_D(const double df_x, const double df_y) ****************************************************************************/ #define NTERM2 10 +//void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax, +// double energy, bool diffusion, bool matlab) void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax, - double energy, bool diffusion, bool matlab) + double energy, bool diffusion) { FILE * outf; const char fic[] = "fmap.out"; @@ -810,10 +832,13 @@ void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax, exit_(1); } - fprintf(outf,"# TRACY II v. 2.6 -- %s -- %s \n", fic, asctime2(newtime)); + fprintf(outf,"# TRACY III SYNCHROTRON SOLEIL-- %s -- %s \n", fic, asctime2(newtime)); fprintf(outf,"# nu = f(x) \n"); - fprintf(outf,"# x[mm] z[mm] fx fz" - " dfx dfz D=log_10(sqrt(df_x^2+df_y^2))\n"); +// fprintf(outf,"# x[mm] z[mm] fx fz" +// " dfx dfz D=log_10(sqrt(df_x^2+df_y^2))\n"); +// + fprintf(outf,"# x[m] z[m] fx fz" + " dfx dfz\n"); if ((Nbx < 1) || (Nbz < 1)) @@ -831,13 +856,16 @@ void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax, zp = zp0; // Tracking part + NAFF - for (i = -Nbx; i <= Nbx; i++) { - x = x0 + sgn(i)*sqrt((double)abs(i))*xstep; - if (!matlab) fprintf(outf,"\n"); + for (i = 0; i <= Nbx; i++) { + x = x0 + sqrt((double)i)*xstep; +// for (i = -Nbx; i <= Nbx; i++) { + // x = x0 + sgn(i)*sqrt((double)abs(i))*xstep; +// if (!matlab) fprintf(outf,"\n"); fprintf(stdout,"\n"); -// for (j = 0; j<= Nbz; j++) { - for (j = -Nbz; j<= Nbz; j++) { - z = z0 + sgn(j)*sqrt((double)abs(j))*zstep; + for (j = 0; j<= Nbz; j++) { + z = z0 + sqrt((double)j)*zstep; + // for (j = -Nbz; j<= Nbz; j++) { + // z = z0 + sgn(j)*sqrt((double)abs(j))*zstep; // tracking around closed orbit Trac_Simple(x,xp,z,zp,energy,ctau,nturn,Tab,&status); if (status) { // if trajectory is stable @@ -859,17 +887,25 @@ void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax, // printout value if (!diffusion){ - fprintf(outf,"%14.6e %14.6e %14.6e %14.6e\n", - 1e3*x, 1e3*z, nux1, nuz1); +// fprintf(outf,"%14.6e %14.6e %14.6e %14.6e\n", +// 1e3*x, 1e3*z, nux1, nuz1); +// fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n", +// 1e3*x, 1e3*z, nux1, nuz1); + fprintf(outf,"%10.6e %10.6e %10.6e %10.6e\n", + x, z, nux1, nuz1); fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n", - 1e3*x, 1e3*z, nux1, nuz1); + x, z, nux1, nuz1); } else { dfx = nux1 - nux2; dfz = nuz1 - nuz2; - fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n", - 1e3*x, 1e3*z, nux1, nuz1, dfx, dfz, get_D(dfx, dfz)); - fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n", - 1e3*x, 1e3*z, nux1, nuz1, dfx, dfz, get_D(dfx, dfz)); +// fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n", +// 1e3*x, 1e3*z, nux1, nuz1, dfx, dfz, get_D(dfx, dfz)); +// fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n", +// 1e3*x, 1e3*z, nux1, nuz1, dfx, dfz, get_D(dfx, dfz)); + fprintf(outf,"%10.6e %10.6e %10.6e %10.6e %10.6e %10.6e\n", + x, z, nux1, nuz1, dfx, dfz); + fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n", + x, z, nux1, nuz1, dfx, dfz); } } } @@ -886,19 +922,23 @@ void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax, Compute a frequency map of Nbx x Nbz points For each set of initial conditions the particle is tracked over Nbtour for an energy offset dp - + + Frequency map is based on fixed vertical amplitude z, trace x versus energy, + or, tracking x for off-momentum particle. + The stepsize follows a square root law Results in fmapdp.out Input: Nbx horizontal step number + Nbe energy step number Nbdp energy step number xmax horizontal maximum amplitude emax energy z vertical starting amplitude Nbtour number of turn for tracking - + matlab set file print format for matlab post-process; specific for nsrl-ii Output: status true if stable false otherwise @@ -919,8 +959,10 @@ void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax, ****************************************************************************/ #define NTERM2 10 +//void fmapdp(long Nbx, long Nbe, long Nbtour, double xmax, double emax, +// double z, bool diffusion, bool matlab) void fmapdp(long Nbx, long Nbe, long Nbtour, double xmax, double emax, - double z, bool diffusion, bool matlab) + double z, bool diffusion) { FILE * outf; const char fic[] = "fmapdp.out"; @@ -952,10 +994,11 @@ void fmapdp(long Nbx, long Nbe, long Nbtour, double xmax, double emax, exit_(1); } - fprintf(outf,"# TRACY II v. 2.6 -- %s -- %s \n", fic, asctime2(newtime)); + fprintf(outf,"# TRACY III SYNCHROTRON SOLEIL-- %s -- %s \n", fic, asctime2(newtime)); fprintf(outf,"# nu = f(x) \n"); - fprintf(outf,"# dp[%%] x[mm] fx fz dfx dfz\n"); - +// fprintf(outf,"# dp[%%] x[mm] fx fz dfx dfz\n"); +fprintf(outf,"# dp[m] x[m] fx fz dfx dfz\n"); + if ((Nbx <= 1) || (Nbe <= 1)) fprintf(stdout,"fmap: Error Nbx=%ld Nbe=%ld\n",Nbx,Nbe); @@ -967,19 +1010,20 @@ void fmapdp(long Nbx, long Nbe, long Nbtour, double xmax, double emax, for (i = 0; i <= Nbe; i++) { dp = -emax + i*estep; - if (!matlab) fprintf(outf,"\n"); +// if (!matlab) fprintf(outf,"\n"); fprintf(stdout,"\n"); -// for (j = 0; j<= Nbx; j++) { - for (j = -Nbx; j<= Nbx; j++) { + for (j = 0; j<= Nbx; j++) { +// for (j = -Nbx; j<= Nbx; j++) { // IF 6D Tracking diffusion turn off and x negative for dp negative if ((globval.Cavity_on == true) && (dp < 0.0)){ - x = x0 - sgn(j)*sqrt((double)abs(j))*xstep; - diffusion = false; + // x = x0 - sgn(j)*sqrt((double)abs(j))*xstep; + x = x0 - sqrt((double)j)*xstep; + diffusion = false; } else - x = x0 + sgn(j)*sqrt((double)abs(j))*xstep; - + // x = x0 + sgn(j)*sqrt((double)abs(j))*xstep; + x = x0 + sqrt((double)j)*xstep; Trac_Simple(x,xp,z,zp,dp,ctau,nturn,Tab,&status); if (status) { Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq); @@ -997,17 +1041,23 @@ void fmapdp(long Nbx, long Nbe, long Nbtour, double xmax, double emax, // printout value if (!diffusion){ - fprintf(outf,"%14.6e %14.6e %14.6e %14.6e\n", - 1e2*dp, 1e3*x, nux1, nuz1); - fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n", - 1e2*dp, 1e3*x, nux1, nuz1); +// fprintf(outf,"%14.6e %14.6e %14.6e %14.6e\n", +// 1e2*dp, 1e3*x, nux1, nuz1); +// fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n", +// 1e2*dp, 1e3*x, nux1, nuz1); + fprintf(outf,"% 10.6e % 10.6e % 10.6e % 10.6e\n", dp, x, nux1, nuz1); + fprintf(stdout,"% 10.6e % 10.6e % 10.6e % 10.6e\n", dp, x, nux1, nuz1); } else { dfx = nux2 - nux1; dfz = nuz2 - nuz1; - fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n", - 1e2*dp, 1e3*x, nux1, nuz2, dfx, dfz, get_D(dfx, dfz)); - fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n", - 1e2*dp, 1e3*x, nux1, nuz2, dfx, dfz, get_D(dfx, dfz)); +// fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n", +// 1e2*dp, 1e3*x, nux1, nuz2, dfx, dfz, get_D(dfx, dfz)); +// fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n", +// 1e2*dp, 1e3*x, nux1, nuz2, dfx, dfz, get_D(dfx, dfz)); + fprintf(outf,"% 10.6e % 10.6e % 10.6e % 10.6e % 10.6e % 10.6e\n", + dp, x, nux1, nuz2, dfx, dfz); + fprintf(stdout,"% 10.6e % 10.6e % 10.6e % 10.6e % 10.6e % 10.6e\n", + dp, x, nux1, nuz2, dfx, dfz); } } } @@ -1071,7 +1121,7 @@ void NuDp(long Nb, long Nbtour, double emax) exit_(1); } - fprintf(outf,"# TRACY II v. 2.6 -- %s -- %s \n", fic, asctime2(newtime)); + fprintf(outf,"# TRACY III -- %s -- %s \n", fic, asctime2(newtime)); fprintf(outf,"# dP/P fx fz xcod pxcod zcod pzcod\n"); fprintf(stdout,"# dP/P fx fz xcod pxcod zcod pzcod\n"); @@ -1166,7 +1216,7 @@ void Phase(double x,double xp,double y, double yp,double energy, double ctau, lo exit_(1); } - fprintf(outf,"# TRACY II v. 2.6 -- %s -- %s \n", fic, asctime2(newtime)); + fprintf(outf,"# TRACY III -- %s -- %s \n", fic, asctime2(newtime)); fprintf(outf,"# Phase Space \n"); fprintf(outf, "# x xp z zp dp ctau\n"); @@ -1254,7 +1304,7 @@ void PhasePoly(long pos, double x0,double px0, double z0, double pz0, double del exit_(1); } - fprintf(outf,"# TRACY II v. 2.6 -- %s -- %s \n", fic, asctime2(newtime)); + fprintf(outf,"# TRACY III -- %s -- %s \n", fic, asctime2(newtime)); fprintf(outf,"# 6D Phase Space \n"); fprintf(outf, "# num x xp z zp dp ctau\n"); @@ -1333,7 +1383,7 @@ void PhasePortrait(double x0,double px0,double z0, double pz0, double delta0, exit_(1); } - fprintf(outf,"# TRACY II v. 2.6 -- %s \n", asctime2(newtime)); + fprintf(outf,"# TRACY III -- %s \n", asctime2(newtime)); fprintf(outf,"# x xp z zp dp ctau\n#\n"); x = x0; px = px0; @@ -1514,11 +1564,872 @@ void Enveloppe(double x, double px, double y, double py, double dp, double nturn } +/****************************************************************************/ +/* void Multipole_new(void) + + Purpose: + Set multipole in dipoles, quadrupoles, thick sextupoles, skew quadrupole, + horizontal and vertical corrector. + + Input: + none + + Output: + none + + Return: + none + + Global variables: + trace + + Specific functions: + getelem, SetKLpar, GetKpar + + Comments: + Test for short and long quadrupole could be changed using the length + instead of the name. Maybe more portable, in particular if periodicity + is broken + Should be rewritten because list already exists now .. + +****************************************************************************/ + +// void Multipole_new(void) +// { +// int i = 0; +// int ndip = 0, /* Number of dipoles */ +// nquad = 0, /* Number of quadrupoles */ +// nsext = 0, /* Number of sextupoles */ +// nhcorr= 0, /* Number of horizontal correctors */ +// nvcorr= 0, /* Number of vertical correctors */ +// nqcorr= 0; /* Number of skew quadrupoles */ + +// int dlist[500]; /* dipole list */ +// int qlist[500]; /* Quadrupole list */ +// int slist[500]; /* Sextupole list */ +// int hcorrlist[120]; /* horizontal corrector list */ +// int vcorrlist[120]; /* vertical corrector list */ +// int qcorrlist[120]; /* skew quad list */ +// int hcorrlistThick[120]; /* horizontal corrector list */ +// int vcorrlistThick[120]; /* vertical corrector list */ +// int qcorrlistThick[120]; /* skew quad list */ + +// CellType Cell; + +// int ElemFamIndex=0;; +// int mOrder = 0; /* multipole order */ +// double mKL = 0.0 ; /* multipole integrated strength */ +// double corr_strength = 0.0; +// double hcorr[120], vcorr[120], qcorr[120]; +// double hcorrThick[120], vcorrThick[120], qcorrThick[120]; +// double b2 = 0.0, b3 = 0.0; +// double dBoB2 = 0.0, dBoB3 = 0.0, dBoB4 = 0.0, dBoB5 = 0.0, dBoB6 = 0.0, +// dBoB7 = 0.0, dBoB9 = 0.0, dBoB11 = 0.0, dBoB15 = 0.0, dBoB21 = 0.0, +// dBoB27; +// double dBoB6C = 0.0, dBoB6L = 0.0, dBoB10C = 0.0, dBoB10L = 0.0, +// dBoB14C = 0.0, dBoB14L = 0.0, dBoB3C = 0.0, dBoB3L = 0.0, +// dBoB4C = 0.0, dBoB4L = 0.0; +// double dBoB5rms = 0.0, dBoB7rms = 0.0; +// double x0i = 0.0, x02i = 0.0, x03i = 0.0, x04i = 0.0, x05i = 0.0, +// x06i = 0.0, x07i = 0.0, x08i = 0.0, x012i = 0.0, x010i = 0.0, +// x018i = 0.0, x024i = 0.0, x1i = 0.0; +// double theta = 0.0, brho = 0.0, dummyf = 0.0, conv = 0.0 ; +// char *dummy; + +// FILE *fi; +// char *fic_hcorr="//home/nadolski/codes/tracy/maille/soleil/corh.txt"; +// char *fic_vcorr="/home/nadolski/codes/tracy/maille/soleil/corv.txt"; +// char *fic_skew ="/home/nadolski/codes/tracy/maille/soleil/corqt.txt"; +// /*********************************************************/ + +// getglobv_(&globval); + +// printf("Enter multipole ... \n"); + +// /* Make lists of dipoles, quadrupoles and sextupoles */ +// for (i = 0; i <= globval.Cell_nLoc; i++) +// { +// getelem(i, &Cell); /* get element */ + +// if (Cell.Elem.Pkind == Mpole) +// { +// if (Cell.Elem.M->UU.U0.Pirho!= 0.0) +// { +// dlist[ndip] = i; +// ndip++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PB[0 + HOMmax]); +// }/* +// else if (Cell.Elem.M->PBpar[2L + HOMmax] != 0.0) +// { +// qlist[nquad] = i; +// nquad++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]); +// } +// else if (Cell.Elem.M->PBpar[3L + HOMmax] != 0.0) +// { +// slist[nsext] = i; +// nsext++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[3L + HOMmax]); +// } */ +// else if ( Cell.Elem.PName[0] == 'c' && Cell.Elem.PName[1] == 'h') +// { +// hcorrlist[nhcorr] = i; +// nhcorr++; +// if (trace) printf("%s \n",Cell.Elem.PName); +// } +// else if ( Cell.Elem.PName[0] == 'c' && Cell.Elem.PName[1] == 'v') +// { +// vcorrlist[nvcorr] = i; +// nvcorr++; +// if (trace) printf("%s \n",Cell.Elem.PName); +// } +// else if ( Cell.Elem.PName[0] == 'q' && Cell.Elem.PName[1] == 't') +// { +// qcorrlist[nqcorr] = i; +// nqcorr++; +// if (trace) printf("%s \n",Cell.Elem.PName); +// } +// } +// } + +// /* building list of quadrupoles */ +// nquad = 0; +// ElemFamIndex = ElemIndex("QP1"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// qlist[nquad] = Elem_GetPos(ElemFamIndex, i); +// getelem(qlist[nquad], &Cell); /* get element */ +// nquad++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("QP2"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// qlist[nquad] = Elem_GetPos(ElemFamIndex, i); +// getelem(qlist[nquad], &Cell); /* get element */ +// nquad++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("QP3"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// qlist[nquad] = Elem_GetPos(ElemFamIndex, i); +// getelem(qlist[nquad], &Cell); /* get element */ +// nquad++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("QP4"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// qlist[nquad] = Elem_GetPos(ElemFamIndex, i); +// getelem(qlist[nquad], &Cell); /* get element */ +// nquad++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("QP5"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// qlist[nquad] = Elem_GetPos(ElemFamIndex, i); +// getelem(qlist[nquad], &Cell); /* get element */ +// nquad++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("QP6"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// qlist[nquad] = Elem_GetPos(ElemFamIndex, i); +// getelem(qlist[nquad], &Cell); /* get element */ +// nquad++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("QP7"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// qlist[nquad] = Elem_GetPos(ElemFamIndex, i); +// getelem(qlist[nquad], &Cell); /* get element */ +// nquad++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("QP8"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// qlist[nquad] = Elem_GetPos(ElemFamIndex, i); +// getelem(qlist[nquad], &Cell); /* get element */ +// nquad++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("QP9"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// qlist[nquad] = Elem_GetPos(ElemFamIndex, i); +// getelem(qlist[nquad], &Cell); /* get element */ +// nquad++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("QP10"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// qlist[nquad] = Elem_GetPos(ElemFamIndex, i); +// getelem(qlist[nquad], &Cell); /* get element */ +// nquad++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[2L + HOMmax]); +// } + +// /* building list of sextupoles */ +// nsext = 0; +// ElemFamIndex = ElemIndex("SX1"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// slist[nsext] = Elem_GetPos(ElemFamIndex, i); +// getelem(slist[nsext], &Cell); /* get element */ +// nsext++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[3L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("SX2"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// slist[nsext] = Elem_GetPos(ElemFamIndex, i); +// getelem(slist[nsext], &Cell); /* get element */ +// nsext++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[3L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("SX3"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// slist[nsext] = Elem_GetPos(ElemFamIndex, i); +// getelem(slist[nsext], &Cell); /* get element */ +// nsext++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[3L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("SX4"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// slist[nsext] = Elem_GetPos(ElemFamIndex, i); +// getelem(slist[nsext], &Cell); /* get element */ +// nsext++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[3L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("SX5"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// slist[nsext] = Elem_GetPos(ElemFamIndex, i); +// getelem(slist[nsext], &Cell); /* get element */ +// nsext++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[3L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("SX6"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// slist[nsext] = Elem_GetPos(ElemFamIndex, i); +// getelem(slist[nsext], &Cell); /* get element */ +// nsext++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[3L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("SX7"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// slist[nsext] = Elem_GetPos(ElemFamIndex, i); +// getelem(slist[nsext], &Cell); /* get element */ +// nsext++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[3L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("SX8"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// slist[nsext] = Elem_GetPos(ElemFamIndex, i); +// getelem(slist[nsext], &Cell); /* get element */ +// nsext++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[3L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("SX9"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// slist[nsext] = Elem_GetPos(ElemFamIndex, i); +// getelem(slist[nsext], &Cell); /* get element */ +// nsext++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[3L + HOMmax]); +// } + +// ElemFamIndex = ElemIndex("SX10"); +// for (i=1; i<=GetnKid(ElemFamIndex); i++){ +// slist[nsext] = Elem_GetPos(ElemFamIndex, i); +// getelem(slist[nsext], &Cell); /* get element */ +// nsext++; +// if (trace) printf("%s % f\n",Cell.Elem.PName, Cell.Elem.M->PBpar[3L + HOMmax]); +// } + + +// /* find sextupole associated with the corrector */ +// // solution 1: find by names +// // solution 2: use a predfined list +// // solution 3: smothing smart ??? +// for (i=0; i< nhcorr; i++){ +// if (trace) fprintf(stdout, "%d\n", i); +// getelem(hcorrlist[i]-1, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// hcorrlistThick[i] = hcorrlist[i]-1; +// else{ +// getelem(hcorrlist[i]+1, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// hcorrlistThick[i] = hcorrlist[i]+1; +// else{ +// getelem(hcorrlist[i]+2, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// hcorrlistThick[i] = hcorrlist[i]+2; +// else{ +// getelem(hcorrlist[i]-2, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// hcorrlistThick[i] = hcorrlist[i]-2; +// else{ +// getelem(hcorrlist[i]+3, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// hcorrlistThick[i] = hcorrlist[i]+3; +// else{ +// getelem(hcorrlist[i]-3, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// hcorrlistThick[i] = hcorrlist[i]-3; +// else fprintf(stdout, "Warning Sextupole not found for VCOR\n"); +// } +// } +// } +// } +// } +// } + +// for (i=0; i< nvcorr; i++){ +// if (trace) fprintf(stdout, "%d\n", i); +// getelem(vcorrlist[i]-1, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// vcorrlistThick[i] = vcorrlist[i]-1; +// else{ +// getelem(vcorrlist[i]+1, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// vcorrlistThick[i] = vcorrlist[i]+1; +// else{ +// getelem(vcorrlist[i]+2, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// vcorrlistThick[i] = vcorrlist[i]+2; +// else{ +// getelem(vcorrlist[i]-2, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// vcorrlistThick[i] = vcorrlist[i]-2; +// else{ +// getelem(vcorrlist[i]+3, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// vcorrlistThick[i] = vcorrlist[i]+3; +// else{ +// getelem(vcorrlist[i]-3, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// vcorrlistThick[i] = vcorrlist[i]-3; +// else fprintf(stdout, "Warning Sextupole not found for VCOR\n"); +// } +// } +// } +// } +// } +// } + +// for (i=0; i< nqcorr; i++){ +// if (trace) fprintf(stdout, "%d\n", i); +// getelem(qcorrlist[i]-1, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// qcorrlistThick[i] = qcorrlist[i]-1; +// else{ +// getelem(qcorrlist[i]+1, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// qcorrlistThick[i] = qcorrlist[i]+1; +// else{ +// getelem(qcorrlist[i]+2, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// qcorrlistThick[i] = qcorrlist[i]+2; +// else{ +// getelem(qcorrlist[i]-2, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// qcorrlistThick[i] = qcorrlist[i]-2; +// else{ +// getelem(qcorrlist[i]+3, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// qcorrlistThick[i] = qcorrlist[i]+3; +// else{ +// getelem(qcorrlist[i]-3, &Cell); +// if (Cell.Elem.PName[0] == 's' && Cell.Elem.PName[1] == 'x') +// qcorrlistThick[i] = qcorrlist[i]-3; +// else fprintf(stdout, "Warning Sextupole not found for QT\n"); +// } +// } +// } +// } +// } +// } + + +// if (!trace) printf("Elements: ndip=%d nquad=%d nsext=%d nhcorr=%d nvcorr=%d nqcorr=%d\n", +// ndip,nquad,nsext,nhcorr,nvcorr,nqcorr); + +// /***********************************************************************************/ +// /* +// /*********** Set multipoles for dipole ****************/ +// /* +// * x0ni w/ n = p-1 for a 2p-poles +// * +// /***********************************************************************************/ + +// x0i = 1.0/20e-3; /* 1/radius */ +// x02i = x0i*x0i; +// x03i = x02i*x0i; +// x04i = x02i*x02i; +// x05i = x04i*x0i; +// x06i = x03i*x03i; +// x07i = x06i*x0i; + +// dBoB2 = 2.2e-4*0; /* gradient, used for curve trajectory simulation */ +// dBoB3 = -3.0e-4*1; /* hexapole */ +// dBoB4 = 2.0e-5*1; /* octupole */ +// dBoB5 = -1.0e-4*1; /* decapole */ +// dBoB6 = -6.0e-5*1; /* 12-poles */ +// dBoB7 = -1.0e-4*1; /* 14-poles */ + +// for (i = 0; i < ndip; i++) +// { +// getelem(dlist[i], &Cell); +// theta = Cell.Elem.PL*Cell.Elem.M->UU.U0.Pirho; + +// /* gradient error */ +// mKL =GetKLpar(Cell.Fnum, Cell.Knum, mOrder=2L); +// mKL += dBoB2*theta*x0i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=2L, mKL); + +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld theta=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, theta, mKL); + +// /* sextupole error */ +// mKL = dBoB3*theta*x02i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=3L, mKL); + +// /* octupole error */ +// mKL = dBoB4*theta*x03i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=4L, mKL); + +// /* decapole error */ +// mKL = dBoB5*theta*x04i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=5L, mKL); + +// /* 12-pole error */ +// mKL = dBoB6*theta*x05i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=6L, mKL); + +// /* 14-pole error */ +// mKL = dBoB7*theta*x06i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=7L, mKL); +// } + +// /***********************************************************************************/ +// /* +// /*********** Set multipoles for quadripole ****************/ +// /* +// * x0ni w/ n = p-2 for a 2p-poles +// * +// /***********************************************************************************/ + +// x0i = 1.0/30e-3; /* 1/Radius in meters */ +// b2 = 0.0; /* Quadrupole strength */ +// x02i = x0i*x0i; +// x04i = x02i*x02i; /* 10-poles */ +// x08i = x04i*x04i; /* 20-poles */ +// x012i= x08i*x04i; /* 28-poles */ + +// dBoB6C = 2.4e-4*1; +// dBoB10C = 0.7e-4*1; +// dBoB14C = 0.9e-4*1; +// dBoB6L = 0.7e-4*1; +// dBoB10L = 1.9e-4*1; +// dBoB14L = 1.0e-4*1; + + +// x1i = 1.0/30e-3; /* rayon reference = 30 mm pour mesure sextupole et octupole*/ +// dBoB3L = 2.9e-4*1; /* sextupole qpole long */ +// dBoB4L = -8.6e-4*1; /* octupole qpole long */ +// dBoB3C = -1.6e-4*1; /* sextupole qpole court */ +// dBoB4C = -3.4e-4*1; /* octupole qpole court */ + + +// for (i = 0; i < nquad; i++) +// { +// getelem(qlist[i], &Cell); +// // b2 = Cell.Elem.PL*GetKpar(Cell.Fnum, Cell.Knum, 2L); +// b2 = GetKLpar(Cell.Fnum, Cell.Knum, 2L); + +// /* 12-pole multipole error */ +// if ((strncmp(Cell.Elem.PName,"qp2",3)==0) || (strncmp(Cell.Elem.PName,"qp7",3)==0)) +// mKL= b2*dBoB6L*x04i; +// else +// mKL= b2*dBoB6C*x04i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=6L, mKL); +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL); + +// /* 20-pole multipole error */ +// if ((strncmp(Cell.Elem.PName,"qp2",3)==0) || (strncmp(Cell.Elem.PName,"qp7",3)==0)) +// mKL= b2*dBoB10L*x08i; +// else +// mKL= b2*dBoB10C*x08i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=10L, mKL); +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL); + +// /* 28-pole multipole error */ +// if ((strncmp(Cell.Elem.PName,"qp2",3)==0) || (strncmp(Cell.Elem.PName,"qp7",3)==0)) +// mKL= b2*dBoB14L*x012i; +// else +// mKL= b2*dBoB14C*x012i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=14L, mKL); +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL); + +// /* sextupole mesure quadrupoles longs*/ +// if ((strncmp(Cell.Elem.PName,"qp2",3)==0) || (strncmp(Cell.Elem.PName,"qp7",3)==0)) +// mKL= b2*dBoB3L*x1i; +// else +// mKL= b2*dBoB3C*x1i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=3L, mKL); +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL); + +// /* octupole mesure quadrupoles longs*/ +// if ((strncmp(Cell.Elem.PName,"qp2",3)==0) || (strncmp(Cell.Elem.PName,"qp7",3)==0)) +// mKL= b2*dBoB4L*x1i*x1i; +// else +// mKL= b2*dBoB4C*x1i*x1i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=4L, mKL); +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL); +// } + +// /***********************************************************************************/ +// /* +// /*********** Set multipoles for sextupole ****************/ +// /* +// * x0ni w/ n = p-3 for a 2p-poles +// * +// /***********************************************************************************/ + +// b3 = 0.0; +// x0i = 1.0/32e-3; +// x02i = x0i*x0i; +// x04i = x02i*x02i; +// x06i = x04i*x02i; /* 18-poles */ +// x012i = x06i*x06i; /* 30-poles */ +// x018i = x012i*x06i; /* 42-poles */ +// x024i = x012i*x012i; /* 54-poles */ + +// /* multipoles from dipolar unallowed component */ +// dBoB5 = 5.4e-4*1; +// dBoB7 = 3.3e-4*1; +// dBoB5rms = 4.7e-4*1; +// dBoB7rms = 2.1e-4*1; + +// /* allowed multipoles */ +// dBoB9 = -4.7e-4*1; +// dBoB15 = -9.0e-4*1; +// dBoB21 = -20.9e-4*1; +// dBoB27 = 0.8e-4*1; +// /* +// dBoB9 = 3.1e-3*1; +// dBoB15 = 5.0e-4*1; +// dBoB21 = -2.0e-2*1; +// dBoB27 = 1.1e-2*1; +// */ +// for (i = 0; i < nsext; i++) +// { +// getelem(slist[i], &Cell); +// b3 = GetKLpar(Cell.Fnum, Cell.Knum, 3L); + +// /* 10-pole multipole error */ +// mKL= b3*dBoB5*x02i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=5L, mKL); +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, b3, mKL); + +// /* 14-pole multipole error */ +// mKL= b3*dBoB7*x04i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=7L, mKL); +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, b3, mKL); + +// /* 18-pole multipole error */ +// mKL= b3*dBoB9*x06i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=9L, mKL); +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, b3, mKL); + +// /* 30-pole multipole error */ +// mKL= b3*dBoB15*x012i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=15L, mKL); +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, b3, mKL); + +// /* 42-pole multipole error */ +// mKL= b3*dBoB21*x018i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=21L, mKL); +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, b3, mKL); + +// /* 54-pole multipole error */ +// mKL= b3*dBoB27*x024i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=27L, mKL); +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b3=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, b3, mKL); +// } + +// /***********************************************************************************/ +// /* +// /****** Set multipoles for sextupole horizontal correctors ****************/ +// /* +// * x0ni w/ n = p-1 for a 2p-poles +// * +// /***********************************************************************************/ +// x0i = 1.0/35e-3; /* 1/radius */ +// x02i = x0i*x0i; +// x03i = x02i*x0i; +// x04i = x02i*x02i; +// x05i = x04i*x0i; +// x06i = x03i*x03i; +// x010i = x05i*x05i; + +// dBoB5 = 0.430*1; /* decapole */ +// dBoB7 = 0.063*1; /* 14-poles */ +// dBoB11 =-0.037*1; /* 22-poles */ + +// brho = globval.Energy/0.299792458; /* magnetic rigidity */ +// // brho = 1.0; /* magnetic rigidity */ +// conv = 8.14e-4; /*conversion des A en T.m*/ + +// /* open H corrector file */ +// if ((fi = fopen(fic_hcorr,"r")) == NULL) +// { +// fprintf(stderr, "Error while opening file %s \n",fic_hcorr); +// exit(1); +// } + +// for (i = 0; i < nhcorr; i++) +// { +// fscanf(fi,"%le \n", &hcorr[i]); +// } +// fclose(fi); /* close H corrector file */ + +// // for (i = 0; i < nhcorr; i++){ +// // getelem(hcorrlist[i], &Cell); +// // corr_strength = hcorr[i]*conv/brho; +// // +// // /* decapole error */ +// // mKL = dBoB5*corr_strength*x04i; +// // SetKLpar(Cell.Fnum, Cell.Knum, mOrder=5L, mKL); +// // +// // if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i, +// // Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL); +// // /* 14-pole error */ +// // mKL = dBoB7*corr_strength*x06i; +// // SetKLpar(Cell.Fnum, Cell.Knum, mOrder=7L, mKL); +// // +// // if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i, +// // Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL); +// // +// // /* 22-pole error */ +// // mKL = dBoB11*corr_strength*x010i; +// // SetKLpar(Cell.Fnum, Cell.Knum, mOrder=11, mKL); +// // +// // if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i, +// // Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL); +// // } + +// for (i = 0; i < nhcorr; i++){ +// getelem(hcorrlistThick[i], &Cell); +// corr_strength = hcorr[i]*conv/brho; + +// /* gradient error */ +// mKL = GetKLpar(Cell.Fnum, Cell.Knum, mOrder=5L); +// mKL += dBoB5*corr_strength*x04i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=5L, mKL); + +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL); +// /* 14-pole error */ +// mKL = GetKLpar(Cell.Fnum, Cell.Knum, mOrder=7L); +// mKL += dBoB7*corr_strength*x06i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=7L, mKL); + +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL); + +// /* 22-pole error */ +// mKL = GetKLpar(Cell.Fnum, Cell.Knum, mOrder=11L); +// mKL += dBoB11*corr_strength*x010i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=11, mKL); + +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL); +// } + +// /***********************************************************************************/ +// /* +// /****** Set multipoles for vertical correctors ****************/ +// /* +// * x0ni w/ n = p-1 for a 2p-poles +// * +// /***********************************************************************************/ + +// x0i = 1.0/35e-3; /* 1/radius */ +// x02i = x0i*x0i; +// x03i = x02i*x0i; +// x04i = x02i*x02i; +// x05i = x04i*x0i; +// x06i = x03i*x03i; +// x010i = x05i*x05i; + +// dBoB5 = -0.430*1; /* decapole */ +// dBoB7 = 0.063*1; /* 14-poles */ +// dBoB11 = 0.037*1; /* 22-poles */ + +// brho = globval.Energy/0.299792458; /* magnetic rigidity */ +// conv = 4.642e-4; /*conversion des A en T.m*/ + + +// /* open V corrector file */ +// if ((fi = fopen(fic_vcorr,"r")) == NULL) +// { +// fprintf(stderr, "Error while opening file %s \n",fic_vcorr); +// exit(1); +// } + +// for (i = 0; i < nvcorr; i++){ +// fscanf(fi,"%le\n", &vcorr[i]); +// } +// fclose(fi); /* close V corrector file */ + +// // for (i = 0; i < nvcorr; i++) +// // { +// // getelem(vcorrlist[i], &Cell); +// // corr_strength = vcorr[i]*conv/brho; +// // +// // /* skew decapole error */ +// // mKL = dBoB5*corr_strength*x04i; +// // SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-5L, mKL); +// // +// // if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i, +// // Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL); +// // /* skew 14-pole error */ +// // mKL = dBoB7*corr_strength*x06i; +// // SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-7L, mKL); +// // +// // if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i, +// // Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL); +// // +// // /* skew 22-pole error */ +// // mKL = dBoB11*corr_strength*x010i; +// // SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-11L, mKL); +// // +// // if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i, +// // Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL); +// // } + +// for (i = 0; i < nvcorr; i++) +// { +// getelem(vcorrlistThick[i], &Cell); +// corr_strength = vcorr[i]*conv/brho; + +// /* skew decapole error */ +// mKL = GetKLpar(Cell.Fnum, Cell.Knum, mOrder=-5L); +// mKL += dBoB5*corr_strength*x04i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-5L, mKL); + +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL); + +// /* skew 14-pole error */ +// mKL = GetKLpar(Cell.Fnum, Cell.Knum, mOrder=-7L); +// mKL += dBoB7*corr_strength*x06i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-7L, mKL); + +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL); + +// /* skew 22-pole error */ +// mKL = GetKLpar(Cell.Fnum, Cell.Knum, mOrder=-11L); +// mKL += dBoB11*corr_strength*x010i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-11L, mKL); + +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL); +// } +// /***********************************************************************************/ +// /* +// /****** Set multipoles for skew quadripole ****************/ +// /* +// * x0ni w/ n = p-2 for a 2p-poles +// * +// /***********************************************************************************/ + +// /* Set multipoles for skew quad */ +// x0i = 1.0/35e-3; /* 1/radius */ +// x02i = x0i*x0i; + +// dBoB4 = -0.680*0; /* Octupole */ + +// /* open skew quaI (A) * +// 310 +// 450 +// 500 +// 520 +// 540 +// 550 +// 560 +// d file */ + +// // brho = 2.75/0.299792458; /* magnetic rigidity */ +// brho = globval.Energy/0.299792458; /* magnetic rigidity */ +// conv = 93.83e-4; /*conversion des A en T*/ + + +// if ((fi = fopen(fic_skew,"r")) == NULL) +// { +// fprintf(stderr, "Error while opening file %s \n",fic_skew); +// exit(1); +// } + +// for (i = 0; i < nqcorr; i++) +// { +// fscanf(fi,"%le \n", &qcorr[i]); +// } +// fclose(fi); /* close skew quad file */ + +// // for (i = 0; i < nqcorr; i++) +// // { +// // getelem(qcorrlist[i], &Cell); +// // +// // /* skew octupole */ +// // mKL = dBoB4*qcorr[i]*conv/brho*x02i; +// // SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-4L, mKL); +// // +// // if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i, +// // Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL); +// // } +// for (i = 0; i < nqcorr; i++) +// { +// getelem(qcorrlist[i], &Cell); + +// /* skew octupole */ +// mKL = GetKLpar(Cell.Fnum, Cell.Knum, mOrder=-4L); +// mKL += dBoB4*qcorr[i]*conv/brho*x02i; +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-4L, mKL); + +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld BL/brho=% e mKl=% e\n",i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, corr_strength, mKL); +// } +// } + /****************************************************************************/ /* void Multipole(void) Purpose: - Set multipole in dipoles, quadrupoles, sextupoles, skew quadrupole, + Set multipole in dipoles, quadrupoles, thin sextupoles, skew quadrupole, horizontal and vertical corrector. Input: @@ -1543,6 +2454,7 @@ void Enveloppe(double x, double px, double y, double py, double dp, double nturn Should be rewritten because list already exists now .. ****************************************************************************/ + void Multipole(void) { int i = 0; @@ -1570,12 +2482,13 @@ void Multipole(void) double dBoB2 = 0.0, dBoB3 = 0.0, dBoB4 = 0.0, dBoB5 = 0.0, dBoB6 = 0.0, dBoB7 = 0.0, dBoB9 = 0.0, dBoB11 = 0.0, dBoB15 = 0.0, dBoB21 = 0.0, dBoB27; - double dBoB6C = 0.0, dBoB6L = 0.0, dBoB10C = 0.0, dBoB10L = 0.0, - dBoB14C = 0.0, dBoB14L = 0.0; + double dBoB6C = 0.0, dBoB6L = 0.0, dBoB10C = 0.0, dBoB10L = 0.0, + dBoB14C = 0.0, dBoB14L = 0.0, dBoB3C = 0.0, dBoB3L = 0.0, + dBoB4C = 0.0, dBoB4L = 0.0; double x0i = 0.0, x02i = 0.0, x03i = 0.0, x04i = 0.0, x05i = 0.0, x06i = 0.0, x07i = 0.0, x08i = 0.0, x012i = 0.0, x010i = 0.0, - x018i = 0.0, x024i = 0.0; - double theta = 0.0, brho = 0.0, dummyf = 0.0 ; + x018i = 0.0, x024i = 0.0, x1i = 0.0; + double theta = 0.0, brho = 0.0, dummyf = 0.0, conv = 0.0 ; char *dummy = NULL; FILE *fi; @@ -1699,19 +2612,27 @@ void Multipole(void) */ /***********************************************************************************/ - x0i = 1.0/35e-3; /* 1/Radius in meters */ + x0i = 1.0/30e-3; /* 1/Radius in meters */ b2 = 0.0; /* Quadrupole strength */ x02i = x0i*x0i; - x04i = x02i*x02i; /* 12-poles */ + x04i = x02i*x02i; /* 10-poles */ x08i = x04i*x04i; /* 20-poles */ x012i= x08i*x04i; /* 28-poles */ - dBoB6C = 1.6e-5*0; - dBoB10C = -1.7e-4*0; - dBoB6L = 6.5e-5*0; - dBoB10L = 1.0e-4*0; - dBoB14L = 1.0e-4*0; - dBoB14C = 1.0e-4*0; + dBoB6C = 2.4e-4*1; + dBoB10C = 0.7e-4*1; + dBoB14C = 0.9e-4*1; + dBoB6L = 0.7e-4*1; + dBoB10L = 1.9e-4*1; + dBoB14L = 1.0e-4*1; + + + x1i = 1.0/30e-3; /* rayon reference = 30 mm pour mesure sextupole et octupole*/ + dBoB3L = 2.9e-4*1; /* sextupole qpole long */ + dBoB4L = -8.6e-4*1; /* octupole qpole long */ + dBoB3C = -1.6e-4*1; /* sextupole qpole court */ + dBoB4C = -3.4e-4*1; /* octupole qpole court */ + for (i = 0; i < nquad; i++) { @@ -1745,6 +2666,24 @@ void Multipole(void) if (trace) printf("num= %4d name = %s Fnum = %3d, Knum=%3d b2=% e mKl=% e\n",i, Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL); +/* sextupole mesure quadrupoles longs*/ + if ((strncmp(Cell.Elem.PName,"qp2",3)==0) || (strncmp(Cell.Elem.PName,"qp7",3)==0)) + mKL= b2*dBoB3L*x1i; + else + mKL= b2*dBoB3C*x1i; + SetKLpar(Cell.Fnum, Cell.Knum, mOrder=3L, mKL); + if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i, + Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL); + + /* octupole mesure quadrupoles longs*/ + if ((strncmp(Cell.Elem.PName,"qp2",3)==0) || (strncmp(Cell.Elem.PName,"qp7",3)==0)) + mKL= b2*dBoB4L*x1i*x1i; + else + mKL= b2*dBoB4C*x1i*x1i; + SetKLpar(Cell.Fnum, Cell.Knum, mOrder=4L, mKL); + if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld b2=% e mKl=% e\n",i, + Cell.Elem.PName,Cell.Fnum, Cell.Knum, b2, mKL); + } /***********************************************************************************/ @@ -1756,17 +2695,24 @@ void Multipole(void) /***********************************************************************************/ b3 = 0.0; - x0i = 1.0/35e-3; + x0i = 1.0/32e-3; x02i = x0i*x0i; x04i = x02i*x02i; x06i = x04i*x02i; /* 18-poles */ x012i = x06i*x06i; /* 30-poles */ x018i = x012i*x06i; /* 42-poles */ x024i = x012i*x012i; /* 54-poles */ - dBoB9 = 8.1e-3*0; - dBoB15 = -1.1e-3*0; - dBoB21 = -1.1e-3*1; - dBoB27 = -1.1e-3*0; + + dBoB9 = -4.7e-4*1; + dBoB15 = -9.0e-4*1; + dBoB21 = -20.9e-4*1; + dBoB27 = 0.8e-4*1 ; +/* + dBoB9 = 3.1e-3*1; + dBoB15 = 5.0e-4*1; + dBoB21 = -2.0e-2*1; + dBoB27 = 1.1e-2*1; +*/ for (i = 0; i < nsext; i++) { @@ -1813,11 +2759,12 @@ void Multipole(void) x06i = x03i*x03i; x010i = x05i*x05i; - dBoB5 = 0.430*0; /* decapole */ - dBoB7 = 0.067*0; /* 14-poles */ - dBoB11 =-0.017*0; /* 22-poles */ + dBoB5 = 0.430*1; /* decapole */ + dBoB7 = 0.063*1; /* 14-poles */ + dBoB11 =-0.037*1; /* 22-poles */ brho = 2.75/0.299792458; /* magnetic rigidity */ + conv = 8.14e-4; /*conversion des A en T.m*/ /* open H corrector file */ if ((fi = fopen(fic_hcorr,"r")) == NULL) @@ -1835,7 +2782,7 @@ void Multipole(void) for (i = 0; i < nhcorr; i++) { getelem(hcorrlist[i], &Cell); - corr_strength = hcorr[i]/brho; + corr_strength = hcorr[i]*conv/brho; /* gradient error */ mKL = dBoB5*corr_strength*x04i; @@ -1874,11 +2821,12 @@ void Multipole(void) x06i = x03i*x03i; x010i = x05i*x05i; - dBoB5 = -0.430*0; /* decapole */ - dBoB7 = 0.063*0; /* 14-poles */ - dBoB11 = 0.037*0; /* 22-poles */ + dBoB5 = -0.430*1; /* decapole */ + dBoB7 = 0.063*1; /* 14-poles */ + dBoB11 = 0.037*1; /* 22-poles */ brho = 2.75/0.299792458; /* magnetic rigidity */ + conv = 4.642e-4; /*conversion des A en T.m*/ /* open V corrector file */ if ((fi = fopen(fic_vcorr,"r")) == NULL) @@ -1889,14 +2837,15 @@ void Multipole(void) for (i = 0; i < nvcorr; i++) { - fscanf(fi,"%s %le %le %le \n", dummy,&dummyf,&dummyf,&vcorr[i]); - } + // fscanf(fi,"%s %le %le %le \n", dummy,&dummyf,&dummyf,&vcorr[i]); + fscanf(fi,"%le\n", &vcorr[i]); +} fclose(fi); /* close V corrector file */ for (i = 0; i < nvcorr; i++) { getelem(vcorrlist[i], &Cell); - corr_strength = vcorr[i]/brho; + corr_strength = vcorr[i]*conv/brho; /* skew decapole error */ mKL = dBoB5*corr_strength*x04i; @@ -1931,7 +2880,21 @@ void Multipole(void) x0i = 1.0/35e-3; /* 1/radius */ x02i = x0i*x0i; - dBoB4 = -0.680*0; /* Octupole */ + dBoB4 = -0.680*1; /* Octupole */ + + /* open skew quaI (A) * +310 +450 +500 +520 +540 +550 +560 +d file */ + + brho = 2.75/0.299792458; /* magnetic rigidity */ + conv = 93.83e-4; /*conversion des A en T*/ + /* open skew quad file */ if ((fi = fopen(fic_skew,"r")) == NULL) @@ -1951,7 +2914,7 @@ void Multipole(void) getelem(qcorrlist[i], &Cell); /* skew octupole */ - mKL = dBoB4*qcorr[i]*x02i; + mKL = dBoB4*qcorr[i]*conv/brho*x02i; SetKLpar(Cell.Fnum, Cell.Knum, mOrder=-4L, mKL); if (trace) printf("num= %4d name = %s Fnum = %3d, Knum=%3d BL/brho=% e mKl=% e\n",i, @@ -2054,6 +3017,67 @@ void SetSkewQuad(void) } } +/****************************************************************************/ +/* void SetDecapole(void) + + Purpose: + Set decapole in horizontal correctors + + Input: + none + + Output: + none + + Return: + none + + Global variables: + trace + + Specific functions: + GetElem, SetKLpar, GetKpar + + Comments: + none + +****************************************************************************/ +// void SetDecapole(void) +// { +// FILE *fi; +// const char fic_deca[] ="/home/nadolski/soltracy/deca.dat"; +// int i; +// double mKL[56]; /* array for skew quad tilt*/ +// CellType Cell; +// long mOrder=5L; + + +// /* open skew quad file */ +// if ((fi = fopen(fic_deca,"r")) == NULL){ +// fprintf(stderr, "Error while opening file %s \n",fic_deca); +// exit(1); +// } + +// /* read decapole strength */ +// for (i = 0; i < globval.hcorr; i++){ +// fscanf(fi,"%le \n", &mKL[i]); +// } +// fclose(fi); + +// for (i = 0; i < globval.hcorr; i++){ +// if (trace) fprintf(stdout,"%le \n", mKL[i]); + +// getelem(globval.hcorr_list[i], &Cell); +// SetKLpar(Cell.Fnum, Cell.Knum, mOrder, mKL[i]); + +// if (trace) printf("num= %4d name = %s Fnum = %3ld, Knum=%3ld KL=% e, KtiltL=% e\n" +// ,i, +// Cell.Elem.PName,Cell.Fnum, Cell.Knum, +// Cell.Elem.M->PBpar[HOMmax+mOrder], +// Cell.Elem.M->PBpar[HOMmax-mOrder]); +// } +// } + /****************************************************************************/ /* void MomentumAcceptance(long deb, long fin, double ep_min, double ep_max, long nstepp, @@ -2082,6 +3106,7 @@ void SetSkewQuad(void) Output: output file soleil.out : file of results + output file phase.out : file of tracking results during the process Return: none @@ -2118,6 +3143,7 @@ void MomentumAcceptance(long deb, long fin, double ep_min, double ep_max, struct tm *newtime; // for time Vector codvector[Cell_nLocMax]; bool cavityflag, radiationflag; + bool trace=true; x0.zero(); @@ -2132,10 +3158,10 @@ void MomentumAcceptance(long deb, long fin, double ep_min, double ep_max, outf1 = fopen("phase.out", "w"); outf2 = fopen("soleil.out", "w"); - fprintf(outf2,"# TRACY II v. 2.6 -- %s \n", asctime2(newtime)); + fprintf(outf2,"# TRACY III v. SYNCHROTRON SOLEIL -- %s \n", asctime2(newtime)); fprintf(outf2,"# i s dp s_lost name_lost \n#\n"); - fprintf(outf1,"# TRACY II v. 2.6 -- %s \n", asctime2(newtime)); + fprintf(outf1,"# TRACY III v. SYNCHROTRON SOLEIL -- %s \n", asctime2(newtime)); fprintf(outf1,"# i x xp z zp dp ctau\n#\n"); @@ -2257,13 +3283,13 @@ void MomentumAcceptance(long deb, long fin, double ep_min, double ep_max, dp2 = ep_max; } if (trace) printf("i=%4ld pos=%4ld dp=%6.4g\n",i,pos,dp2); - if (0) fprintf(stdout,"pos=%4ld z0 =% 10.5f pz0 =% 10.5f \n", pos, tabz0[i-1L][pos-1L], tabpz0[i-1L][pos-1L]); + if (1) fprintf(stdout,"pos=%4ld z0 =% 10.5f pz0 =% 10.5f \n", pos, tabz0[i-1L][pos-1L], tabpz0[i-1L][pos-1L]); Trac(x, px, tabz0[i-1L][pos], tabpz0[i-1L][pos-1L], dp2+delta , ctau0, nturn, pos, lastn, lastpos, outf1); } while (((lastn) == nturn) && (i != nstepp)); - + if ((lastn) == nturn) dp1 = dp2; - + getelem(lastpos,&Clost); getelem(pos,&Cell); fprintf(stdout,"pos=%4ld z0 =% 10.5f pz0 =% 10.5f \n", pos, tabz0[i-1L][pos-1L], tabpz0[i-1L][pos-1L]); @@ -2406,7 +3432,7 @@ void MomentumAcceptance(long deb, long fin, double ep_min, double ep_max, else { dp2 = em_max; } - if (!trace) printf("i=%4ld pos=%4ld dp=%6.4g\n",i,pos,dp2); + if (trace) printf("i=%4ld pos=%4ld dp=%6.4g\n",i,pos,dp2); Trac(x, px, tabz0[i-1L][pos], tabpz0[i-1L][pos-1L], dp2+delta , ctau0, nturn, pos, lastn, lastpos, outf1); } while (((lastn) == nturn) && (i != nstepm)); diff --git a/tracy/tracy/src/t2elem.cc b/tracy/tracy/src/t2elem.cc index fb70019..bfe5e43 100644 --- a/tracy/tracy/src/t2elem.cc +++ b/tracy/tracy/src/t2elem.cc @@ -19,6 +19,14 @@ double **C_; // Dynamic model +/****************************************************************************/ +/* void GtoL(ss_vect<T> &X, Vector2 &S, Vector2 &R, + const double c0, const double c1, const double s1) + + Purpose: + Global to local coordinates + +****************************************************************************/ template<typename T> void GtoL(ss_vect<T> &X, Vector2 &S, Vector2 &R, const double c0, const double c1, const double s1) @@ -39,7 +47,14 @@ void GtoL(ss_vect<T> &X, Vector2 &S, Vector2 &R, X[px_] -= c0; } +/****************************************************************************/ +/* void LtoG(ss_vect<T> &X, Vector2 &S, Vector2 &R, + double c0, double c1, double s1) + Purpose: + Local to global coordinates + +****************************************************************************/ template<typename T> void LtoG(ss_vect<T> &X, Vector2 &S, Vector2 &R, double c0, double c1, double s1) @@ -81,6 +96,38 @@ inline T get_p_s(ss_vect<T> &x) return(p_s); } +/****************************************************************************/ +/* Drift(double L, ss_vect<T> &x) + + Purpose: + Drift pass + + If H_exact = false, using approximation Hamiltonian: + + px^2+py^2 + H(x,y,l,px,py,delta) = ------------- + 2(1+delta) + + Otherwise, using exact Hamiltonian + + + Input: + L: drift length + &x: pointer to initial vector: x + + Output: + none + + Return: + none + + Global variables: + + + Specific functions: + +****************************************************************************/ + template<typename T> void Drift(double L, ss_vect<T> &x) @@ -102,7 +149,14 @@ void Drift(double L, ss_vect<T> &x) template<typename T> void Drift_Pass(CellType &Cell, ss_vect<T> &x) { Drift(Cell.Elem.PL, x); } +/****************************************************************************/ +/* zero_mat(const int n, double** A) + Purpose: + Initionize n x n matrix A with 0 + + +****************************************************************************/ void zero_mat(const int n, double** A) { int i, j; @@ -112,7 +166,14 @@ void zero_mat(const int n, double** A) A[i][j] = 0.0; } +/****************************************************************************/ +/* void identity_mat(const int n, double** A) + Purpose: + generate n x n identity matrix A + + +****************************************************************************/ void identity_mat(const int n, double** A) { int i, j; @@ -122,7 +183,14 @@ void identity_mat(const int n, double** A) A[i][j] = (i == j)? 1.0 : 0.0; } +/****************************************************************************/ +/* void det_mat(const int n, double **A) + Purpose: + get the determinant of n x n matrix A + + +****************************************************************************/ double det_mat(const int n, double **A) { int i, *indx; @@ -140,7 +208,14 @@ double det_mat(const int n, double **A) return d; } +/****************************************************************************/ +/* void trace_mat(const int n, double **A) + Purpose: + get the trace of n x n matrix A + + +****************************************************************************/ double trace_mat(const int n, double **A) { int i; @@ -494,56 +569,100 @@ static double get_psi(double irho, double phi, double gap) } -template<typename T> -void thin_kick(int Order, double MB[], double L, double h_bend, double h_ref, - ss_vect<T> &x) -{ - /* The kick is given by +/****************************************************************************/ +/* template<typename T> +void thin_kick(int Order, double MB[], double L, double h_bend, double h_ref,ss_vect<T> &x) + + Purpose: + Calculate multipole kick. The Hamiltonian is + + H = A + B where A and B are the kick part defined by + 2 2 + px + py + A(x,y,-l,px,py,dP) = --------- + 2(1+dP) + 2 2 + B(x,y,-l,px,py,dP) = -h*x*dP + 0.5 h x + int(Re(By+iBx)/Brho) + + so this is the appproximation for large ring + the chromatic term is missing hx*A + + + The kick is given by e L L delta L x e L Dp_x = - --- B_y + ------- - ----- , Dp_y = --- B_x p_0 rho rho^2 p_0 where - + e 1 + --- = ----- + p_0 B rho ==== \ (B_y + iB_x) = B rho > (ia_n + b_n ) (x + iy)^n-1 / ==== + + Input: + Order maximum non zero multipole component + MB array of an and bn + h_bend 1/rho + h_ref thick element + x initial coordinates vector - where + Output: + x final coordinates vector - e 1 - --- = ----- - p_0 B rho */ + Return: + none + + Global variables: + none + Specific functions: + none + + Comments: + none + +**************************************************************************/ +template<typename T> +void thin_kick(int Order, double MB[], double L, double h_bend, double h_ref, + ss_vect<T> &x) +{ int j; T BxoBrho, ByoBrho, ByoBrho1, B[3]; ss_vect<T> x0, cod; - if ((h_bend != 0.0) || ((1 <= Order) && (Order <= HOMmax))) { - x0 = x; - /* compute field with Horner's rule */ - ByoBrho = MB[Order+HOMmax]; BxoBrho = MB[HOMmax-Order]; - for (j = Order-1; j >= 1; j--) { - ByoBrho1 = x0[x_]*ByoBrho - x0[y_]*BxoBrho + MB[j+HOMmax]; - BxoBrho = x0[y_]*ByoBrho + x0[x_]*BxoBrho + MB[HOMmax-j]; - ByoBrho = ByoBrho1; - } + if ((h_bend != 0.0) || ((1 <= Order) && (Order <= HOMmax))) + { + x0 = x; + /* compute field with Horner's rule */ + ByoBrho = MB[HOMmax+Order]; + BxoBrho = MB[HOMmax-Order]; + + for (j = Order-1; j >= 1; j--) { + ByoBrho1 = x0[x_]*ByoBrho - x0[y_]*BxoBrho + MB[j+HOMmax]; + BxoBrho = x0[y_]*ByoBrho + x0[x_]*BxoBrho + MB[HOMmax-j]; + ByoBrho = ByoBrho1; + } - if (globval.radiation || globval.emittance) { - B[X_] = BxoBrho; B[Y_] = ByoBrho + h_bend; B[Z_] = 0.0; - radiate(x, L, h_ref, B); - } + if (globval.radiation || globval.emittance) { + B[X_] = BxoBrho; + B[Y_] = ByoBrho + h_bend; + B[Z_] = 0.0; + radiate(x, L, h_ref, B); + } - if (h_ref != 0.0) { - x[px_] -= L*(ByoBrho+(h_bend-h_ref)/2.0+h_ref*h_bend*x0[x_] + if (h_ref != 0.0) { + x[px_] -= L*(ByoBrho+(h_bend-h_ref)/2.0+h_ref*h_bend*x0[x_] -h_ref*x0[delta_]); - x[ct_] += L*h_ref*x0[x_]; - } else - x[px_] -= L*(ByoBrho+h_bend); - x[py_] += L*BxoBrho; + x[ct_] += L*h_ref*x0[x_]; + } + else x[px_] -= L*(ByoBrho+h_bend); + + x[py_] += L*BxoBrho; } } @@ -551,13 +670,16 @@ void thin_kick(int Order, double MB[], double L, double h_bend, double h_ref, template<typename T> static void EdgeFocus(double irho, double phi, double gap, ss_vect<T> &x) { - x[px_] += irho*tan(dtor(phi))*x[x_]; - if (false) - // warning: => diverging Taylor map (see SSC-141) - x[py_] -= irho*tan(dtor(phi)-get_psi(irho, phi, gap))*x[y_] + if (true){ + // warning: => diverging Taylor map (see SSC-141) + x[px_] += irho*tan(dtor(phi))*x[x_]/(1.0+x[delta_]); + x[py_] -= irho*tan(dtor(phi)-get_psi(irho, phi, gap))*x[y_] /(1.0+x[delta_]); - else - x[py_] -= irho*tan(dtor(phi)-get_psi(irho, phi, gap))*x[y_]; + } + else{ + x[px_] += irho*tan(dtor(phi))*x[x_]; + x[py_] -= irho*tan(dtor(phi)-get_psi(irho, phi, gap))*x[y_]; + } } @@ -593,68 +715,171 @@ void bend_fringe(double hb, ss_vect<T> &x) } } +/****************************************************************************/ +/* template<typename T> + void quad_fringe(double b2, ss_vect<T> &x) + + Purpose: + Compute edge focusing for a quadrupole + There is no radiation coming from the edge + + The standard formula used is using more general form with exponential + form. If keep up to the 4-th order nonlinear terms, the formula with goes to the + following: + (see E. Forest's book (Beam Dynamics: A New Attitude and Framework), p390): + b2 + x = x0 +/- ---------- (x0^3 + 3*z0^2*x0) + 12(1 + dP) + + b2 + px = px0 +/- ---------- (2*x0*y0*py0 - x0^2*px0 - y0^2*py0) + 4(1 + dP) + + b2 + y = y0 -/+ ---------- (y0^3 + 3*x0^2*y0) + 12(1 + dP) + + b2 + py = py0 -/+ ---------- (2*x0*y0*px0 - y0^2*py0 - x0^2*py0) + 4(1 + dP) + + dP = dP0; + + + b2 + tau = tau0 -/+ ----------- (y0^3*px - x0^3*py + 3*x0^2*y*py - 3*y0^2*x0*px) + 12(1 + dP)^2 + + Input: + b2 = quadrupole strength + x = input particle coordinates + + Output: + x output particle coordinates + + Return: + none + + Global variables: + none + Specific functions: + none + + Comments: + Now in Tracy III, no definition "entrance" and "exit", when called in Mpole_pass, + first call with M --> PB[quad+HOMmax], then + call with -M --> PB[quad+HOMmax] + +****************************************************************************/ template<typename T> void quad_fringe(double b2, ss_vect<T> &x) { T u, ps; - u = b2/(12.0*(1.0+x[delta_])); ps = u/(1.0+x[delta_]); - x[py_] /= 1.0 - 3.0*u*sqr(x[y_]); x[y_] -= u*cube(x[y_]); - if (globval.Cavity_on) x[ct_] -= ps*cube(x[y_])*x[py_]; - x[px_] /= 1.0 + 3.0*u*sqr(x[x_]); - if (globval.Cavity_on) x[ct_] += ps*cube(x[x_])*x[px_]; - x[x_] += u*cube(x[x_]); u = u*3.0; ps = ps*3.0; - x[y_] = exp(-u*sqr(x[x_]))*x[y_]; x[py_] = exp(u*sqr(x[x_]))*x[py_]; - x[px_] += 2.0*u*x[x_]*x[y_]*x[py_]; - if (globval.Cavity_on) x[ct_] -= ps*sqr(x[x_])*x[y_]*x[py_]; - x[x_] = exp(u*sqr(x[y_]))*x[x_]; x[px_] = exp(-u*sqr(x[y_]))*x[px_]; - x[py_] -= 2.0*u*x[y_]*x[x_]*x[px_]; - if (globval.Cavity_on) x[ct_] += ps*sqr(x[y_])*x[x_]*x[px_]; + u = b2/(12.0*(1.0+x[delta_])); + ps = u/(1.0+x[delta_]); + + x[py_] /= 1.0 - 3.0*u*sqr(x[y_]); + x[y_] -= u*cube(x[y_]); + + if (globval.Cavity_on) x[ct_] -= ps*cube(x[y_])*x[py_]; //-y^3*py + x[px_] /= 1.0 + 3.0*u*sqr(x[x_]); //+x^2 + + if (globval.Cavity_on) x[ct_] += ps*cube(x[x_])*x[px_]; //+x^3*px + x[x_] += u*cube(x[x_]); //+x^3 + u = u*3.0; ps = ps*3.0; + x[y_] = exp(-u*sqr(x[x_]))*x[y_]; //+x^2*y + x[py_] = exp(u*sqr(x[x_]))*x[py_]; //+x^2*py + x[px_]+= 2.0*u*x[x_]*x[y_]*x[py_]; //+2*x*y*py + + if (globval.Cavity_on) x[ct_] -= ps*sqr(x[x_])*x[y_]*x[py_]; // -3*x^2*y*py + x[x_] = exp(u*sqr(x[y_]))*x[x_]; //+x*y^2 + x[px_] = exp(-u*sqr(x[y_]))*x[px_]; // -x^2*px-y^2*px + x[py_]-= 2.0*u*x[y_]*x[x_]*x[px_]; //-2*x*y*px + + if (globval.Cavity_on) x[ct_] += ps*sqr(x[y_])*x[x_]*x[px_]; // +3*y^2*x*px } +/****************************************************************************/ +/* void Mpole_Pass(CellType &Cell, ss_vect<T> &x) + + Purpose: + multipole pass,for dipole, quadrupole,sextupole,decupole,etc + Using DA method. + + Input: + + Output: + + + Return: + none + + Global variables: + none + + Specific functions: + none + + Comments: + none + +****************************************************************************/ + template<typename T> void Mpole_Pass(CellType &Cell, ss_vect<T> &x) { - int seg = 0; - double k = 0.0, dL = 0.0, dL1 = 0.0, dL2 = 0.0; + int seg = 0; + double k = 0.0, dL = 0.0, dL1 = 0.0, dL2 = 0.0; double dkL1 = 0.0, dkL2 = 0.0, h_ref = 0.0; elemtype *elemp; MpoleType *M; - + elemp = &Cell.Elem; M = elemp->M; - + /* Global -> Local */ GtoL(x, Cell.dS, Cell.dT, M->Pc0, M->Pc1, M->Ps1); - if ((M->Pmethod == Meth_Second) || (M->Pmethod == Meth_Fourth)) { - /* fringe fields */ + if ((M->Pmethod == Meth_Second) || (M->Pmethod == Meth_Fourth)) + { /* fringe fields */ + if (globval.quad_fringe && (M->PB[Quad+HOMmax] != 0.0)) quad_fringe(M->PB[Quad+HOMmax], x); - if (!globval.H_exact) { - if (M->Pirho != 0.0) EdgeFocus(M->Pirho, M->PTx1, M->Pgap, x); - } else { - p_rot(M->PTx1, x); bend_fringe(M->Pirho, x); - } + + if (!globval.H_exact) + { + if (M->Pirho != 0.0) + EdgeFocus(M->Pirho, M->PTx1, M->Pgap, x); + } + else + { + p_rot(M->PTx1, x); + bend_fringe(M->Pirho, x); + } } - if (M->Pthick == thick) { - if (!globval.H_exact || ((M->PTx1 == 0.0) && (M->PTx2 == 0.0))) { - // polar coordinates - h_ref = M->Pirho; dL = elemp->PL/M->PN; - } else { - // Cartesian coordinates + + if (M->Pthick == thick) + { + if (!globval.H_exact || ((M->PTx1 == 0.0) && (M->PTx2 == 0.0))) + {// polar coordinates,curvilinear coordinates + h_ref = M->Pirho; + dL = elemp->PL/M->PN; + } + else + {// Cartesian coordinates h_ref = 0.0; if (M->Pirho == 0.0) dL = elemp->PL/M->PN; else - dL = 1.0/M->Pirho*(sin(dtor(M->PTx1)) - + sin(dtor(M->PTx2)))/M->PN; + dL = 1.0/M->Pirho*(sin(dtor(M->PTx1)) + sin(dtor(M->PTx2)))/M->PN; } } - switch (M->Pmethod) { + switch (M->Pmethod) + { case Meth_Linear: @@ -682,13 +907,17 @@ void Mpole_Pass(CellType &Cell, ss_vect<T> &x) case Meth_Fourth: if (M->Pthick == thick) { - dL1 = c_1*dL; dL2 = c_2*dL; dkL1 = d_1*dL; dkL2 = d_2*dL; + dL1 = c_1*dL; + dL2 = c_2*dL; + dkL1 = d_1*dL; + dkL2 = d_2*dL; - dcurly_H = 0.0; dI4 = 0.0; + dcurly_H = 0.0; + dI4 = 0.0; for (seg = 1; seg <= M->PN; seg++) { if (globval.emittance && (!globval.Cavity_on) && (M->Pirho != 0.0)) { dcurly_H += is_tps<tps>::get_curly_H(x); - dI4 += is_tps<tps>::get_dI4(x); + dI4 += is_tps<tps>::get_dI4(x); } Drift(dL1, x); @@ -698,7 +927,7 @@ void Mpole_Pass(CellType &Cell, ss_vect<T> &x) if (globval.emittance && (!globval.Cavity_on) && (M->Pirho != 0.0)) { dcurly_H += 4.0*is_tps<tps>::get_curly_H(x); - dI4 += 4.0*is_tps<tps>::get_dI4(x); + dI4 += 4.0*is_tps<tps>::get_dI4(x); } Drift(dL2, x); @@ -707,13 +936,13 @@ void Mpole_Pass(CellType &Cell, ss_vect<T> &x) if (globval.emittance && (!globval.Cavity_on) && (M->Pirho != 0.0)) { dcurly_H += is_tps<tps>::get_curly_H(x); - dI4 += is_tps<tps>::get_dI4(x); + dI4 += is_tps<tps>::get_dI4(x); } } if (globval.emittance && (!globval.Cavity_on) && (M->Pirho != 0)) { dcurly_H /= 6.0*M->PN; - dI4 *= M->Pirho*(sqr(M->Pirho)+2.0*M->PBpar[Quad+HOMmax])/(6.0*M->PN); + dI4 *= M->Pirho*(sqr(M->Pirho)+2.0*M->PBpar[Quad+HOMmax])/(6.0*M->PN); I2 += elemp->PL*sqr(M->Pirho); I4 += elemp->PL*dI4; I5 += elemp->PL*fabs(cube(M->Pirho))*dcurly_H; @@ -724,11 +953,17 @@ void Mpole_Pass(CellType &Cell, ss_vect<T> &x) break; } - if ((M->Pmethod == Meth_Second) || (M->Pmethod == Meth_Fourth)) { + + if ((M->Pmethod == Meth_Second) || (M->Pmethod == Meth_Fourth)) + { /* fringe fields */ - if (!globval.H_exact) { - if (M->Pirho != 0.0) EdgeFocus(M->Pirho, M->PTx2, M->Pgap, x); - } else { + if (!globval.H_exact) + { + if (M->Pirho != 0.0) + EdgeFocus(M->Pirho, M->PTx2, M->Pgap, x); + } + else + { bend_fringe(-M->Pirho, x); p_rot(M->PTx2, x); } if (globval.quad_fringe && (M->PB[Quad+HOMmax] != 0.0)) @@ -2137,8 +2372,10 @@ void SI_init(void) double C_gamma, C_u; - c_1 = 1e0/(2e0*(2e0-thirdroot(2e0))); c_2 = 0.5e0 - c_1; - d_1 = 2e0*c_1; d_2 = 1e0 - 2e0*d_1; + c_1 = 1e0/(2e0*(2e0-thirdroot(2e0))); + c_2 = 0.5e0 - c_1; + d_1 = 2e0*c_1; + d_2 = 1e0 - 2e0*d_1; // classical radiation // C_gamma = 8.846056192e-05; diff --git a/tracy/tracy/src/t2lat.cc b/tracy/tracy/src/t2lat.cc index 896ad82..c948efe 100644 --- a/tracy/tracy/src/t2lat.cc +++ b/tracy/tracy/src/t2lat.cc @@ -24,9 +24,11 @@ J. Bengtsson NSLS-II, BNL 2004 - #define nmax LONG_MAX -#define nn 3 -#define tmax 100 +#define nn 3 // added by nsrl-ii +#define tmax 100 // added by nsrl-ii +#define smax 600 /*size of string table*/ +#define xmax LONG_MAX /*2**8 - 1*/ typedef char Latlinetype[LatLLng]; @@ -4134,10 +4136,16 @@ bool Lattice_Read(FILE **fi_, FILE **fo_) V.fo = fo_; /* output lattice file */ if (setjmp(V._JL9999)) goto _L9999; - V.UDIC = 0; globval.Cell_nLoc = 0; globval.Elem_nFam = 0; V.NoB = 0; - V.Symmetry = 0; V.Bpointer = 0; + V.UDIC = 0; + globval.Cell_nLoc = 0; + globval.Elem_nFam = 0; + V.NoB = 0; + V.Symmetry = 0; + V.Bpointer = 0; - globval.CODeps = 0.0; globval.dPcommon = 0.0; globval.Energy = 0.0; + globval.CODeps = 0.0; + globval.dPcommon = 0.0; + globval.Energy = 0.0; ErrFlag = false; @@ -4145,9 +4153,11 @@ bool Lattice_Read(FILE **fi_, FILE **fo_) GetSym___(&V); - if (V.sym == defsym) DealWithDefns(&V); + if (V.sym == defsym) + DealWithDefns(&V); - if (V.Symmetry != 0) { + if (V.Symmetry != 0) + { GetRingType(&V);/* define whether a ring or a transfer line */ GetEnergy(&V); /* define particle energy */ GetCODEPS(&V); /* define COD precision */ @@ -4176,5 +4186,9 @@ bool Lattice_Read(FILE **fi_, FILE **fo_) #undef kmax #undef nmax -#undef nn -#undef tmax +#undef nn // added by nsrl-ii +#undef tmax // added by nsrl-ii + +#undef smax +#undef xmax + diff --git a/tracy/tracy/src/t2ring.cc b/tracy/tracy/src/t2ring.cc index a00a791..5267f0a 100644 --- a/tracy/tracy/src/t2ring.cc +++ b/tracy/tracy/src/t2ring.cc @@ -391,8 +391,10 @@ void Cell_Twiss(long i0, long i1, ss_vect<tps> &Ascr, bool chroma, bool ring, CellType *cellp; /* initialization */ - for (j = 0; j <= 1; j++) { - nu1[j] = 0.0; dnu[j] = 0.0; + for (j = 0; j <= 1; j++) + { + nu1[j] = 0.0; + dnu[j] = 0.0; } if (globval.radiation) globval.dE = 0.0; @@ -406,8 +408,10 @@ void Cell_Twiss(long i0, long i1, ss_vect<tps> &Ascr, bool chroma, bool ring, Ascr0[j] += tps(globval.CODvect[j]); Ascr1 = Ascr0; - for (i = i0; i <= i1; i++) { - Elem_Pass(i, Ascr1); cellp = &Cell[i]; + for (i = i0; i <= i1; i++) + { + Elem_Pass(i, Ascr1); + cellp = &Cell[i]; dagetprm(Ascr1, cellp->Alpha, cellp->Beta); for (j = 0; j <= 1; j++) { k = (j+1)*2 - 1; @@ -593,7 +597,8 @@ void Ring_Twiss_M(bool chroma, double dP) LinTrans((long)n, C, eta0); for (j = 0; j <= n; j++) { - globval.Ascr[n][j] = 0.0; globval.Ascr[j][n] = eta0[j]; + globval.Ascr[n][j] = 0.0; + globval.Ascr[j][n] = eta0[j]; } CopyMat(n+1, globval.Ascr, Ascr); @@ -662,8 +667,11 @@ void Ring_Twiss(bool chroma, double dP) globval.OneTurnMat, globval.Omega, globval.Alphac); putlinmat(n, globval.Ascr, AScr); - if (!globval.Cavity_on) { - AScr[delta_] = 0.0; AScr[ct_] = 0.0; + + if (!globval.Cavity_on) + { + AScr[delta_] = 0.0; + AScr[ct_] = 0.0; } Cell_Twiss(0, globval.Cell_nLoc, AScr, chroma, true, dP); @@ -682,7 +690,7 @@ void Ring_Twiss(bool chroma, double dP) Purpose: Computes Twiss functions w/ or w/o chromaticities for particle of energy offset dP - matrix or da method + matrix or DA method Input: Chroma if true compute chromaticities and dispersion diff --git a/tracy/tracy/src/tracy.cc b/tracy/tracy/src/tracy.cc index bc86987..f8d82a5 100644 --- a/tracy/tracy/src/tracy.cc +++ b/tracy/tracy/src/tracy.cc @@ -51,6 +51,7 @@ #include "nsls-ii_lib.cc" #include "max4_lib.cc" +#include "soleilcommon.cc" // Truncated Power Series Algebra (TPSA) const int nv_tps = ss_dim, // no of variables -- GitLab