From e2bfe29dc8f723187b5551a7eb19abd97d0eac80 Mon Sep 17 00:00:00 2001
From: nadolski <nadolski@9a6e40ed-f3a0-4838-9b4a-bf418f78e88d>
Date: Fri, 1 Oct 2010 12:53:13 +0000
Subject: [PATCH] removed findcod (was used in Tracy II since da part of
insertion was badly coded)
---
tracy/tracy/src/naffutils.cc | 6 +-
tracy/tracy/src/physlib.cc | 2 +-
tracy/tracy/src/soleilcommon.cc | 809 --------------------------------
tracy/tracy/src/soleillib.cc | 31 +-
4 files changed, 19 insertions(+), 829 deletions(-)
diff --git a/tracy/tracy/src/naffutils.cc b/tracy/tracy/src/naffutils.cc
index db8a2e6..3a6fac0 100644
--- a/tracy/tracy/src/naffutils.cc
+++ b/tracy/tracy/src/naffutils.cc
@@ -45,7 +45,6 @@ double pi = M_PI;
Comments:
useful for connection with NAFF
19/01/03 tracking around the closed orbit
- 20/03/04 closed orbit searched w/ findcod (Newton-Raphson numerical method)
****************************************************************************/
void Trac_Simple(double x, double px, double y, double py, double dp,
double ctau, long nmax, double Tx[][NTURN], bool *status2)
@@ -61,8 +60,9 @@ void Trac_Simple(double x, double px, double y, double py, double dp,
if (globval.MatMeth) Cell_Concat(dp);
/* Get closed orbit */
- //~ getcod(dp, lastpos);
- findcod(dp);
+
+ getcod(dp, lastpos);
+
if (trace && status.codflag)
printf("dp= % .5e %% xcod= % .5e mm zcod= % .5e mm \n",
diff --git a/tracy/tracy/src/physlib.cc b/tracy/tracy/src/physlib.cc
index 3941d70..d3f392f 100644
--- a/tracy/tracy/src/physlib.cc
+++ b/tracy/tracy/src/physlib.cc
@@ -3291,7 +3291,7 @@ void Newton_RaphsonS(int ntrial, double x[], int n, double tolx)
break;
}
}
- // check whever closed orbit found out
+ // check whenver 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);
diff --git a/tracy/tracy/src/soleilcommon.cc b/tracy/tracy/src/soleilcommon.cc
index 31ccdf1..8f8e53e 100644
--- a/tracy/tracy/src/soleilcommon.cc
+++ b/tracy/tracy/src/soleilcommon.cc
@@ -681,816 +681,7 @@ long get_qt_number(void)
return nb;
}
-/****************************************************************************/
-/* 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:
- none
-
-****************************************************************************/
-// void TransTwiss(double *alpha, double *beta, double *eta, double *etap,
-// double *codvect)
-// {
-// TransTrace(0L, globval.Cell_nLoc, alpha, beta, eta, etap, codvect);
-// }
-
-
-
-// /****************************************************************************/
-// /* 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)
-// {
-// float *vcod, *vector();
-// double x0[6];
-// const int ntrial = 40; // maximum number of trials for closed orbit
-// const float tolx = 1e-8; // numerical precision
-// int k;
-// int dim; // 4D or 6D tracking
-// long lastpos;
-// void free_vector();
-// CellType Cell;
-
-// getelem(1,&Cell);
-
-// vcod = vector(1, 6);
-
-// // starting point
-// for (k = 1; k <= 6; k++)
-// vcod[k] = (float) 0.0;
-
-// vcod[5] = (float) dP; // energy offset
-
-// if (globval.Cavity_on){
-// dim = 6; /* 6D tracking*/
-// fprintf(stderr,"Error looking for cod in 6D\n");
-// exit(1);
-// }
-// else{
-// dim = 4; /* 4D tracking */
-// vcod[1] = (float) Cell.Eta[0]*dP;
-// vcod[2] = (float) Cell.Etap[0]*dP;
-// vcod[3] = (float) Cell.Eta[1]*dP;
-// vcod[4] = (float) Cell.Etap[1]*dP;
-// }
-
-// Newton_RaphsonS(ntrial, vcod, dim, (float) tolx);
-
-// if (status.codflag == false)
-// fprintf(stdout, "Error No COD found\n");
-// if (trace)
-// {
-// for (k = 1; k <= 6; k++)
-// x0[k - 1] = (double) 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(1L, 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(1L, globval.Cell_nLoc, x0, &lastpos);
-// }
-// free_vector(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)
-// {
-// double vcod[6];
-// const int ntrial = 40; // maximum number of trials for closed orbit
-// const double tolx = 1e-10; // numerical precision
-// int k, dim = 0;
-// long lastpos;
-// CellType Cell;
-
-// getelem(1,&Cell);
-
-// // initializations
-// for (k = 0; k <= 5; k++)
-// vcod[k] = 0.0;
-
-// if (globval.Cavity_on){
-// fprintf(stderr,"warning looking for cod in 6D\n");
-// dim = 6;
-// }
-// else{ // starting point linear closed orbit
-// dim = 4;
-// vcod[0] = Cell.Eta[0]*dP;
-// vcod[1] = Cell.Etap[0]*dP;
-// vcod[2] = Cell.Eta[1]*dP;
-// vcod[3] = Cell.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(6L, 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(1L, 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(float *x, int n, float **fjac, float *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(float *x, int n, float **fjac, float *fvect)
-// {
-// int i, k;
-// long lastpos = 0L;
-// double x0[6], fx[6], fx1[6], fx2[6];
-// const double deps = 1e-8; //stepsize for numerical differentiation
-
-// for (i = 1; i <= 6; i++)
-// x0[i - 1] = (double) x[i];
-
-// Cell_Pass(1L, globval.Cell_nLoc, x0, &lastpos);
-
-// for (i = 1; i <= n; i++)
-// {
-// fvect[i] = (float) 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] = (double) x[i];
-// x0[k] += deps; // differential step in coordinate k
-
-// Cell_Pass(1L, 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] = (double) x[i];
-// x0[5] = 0.0;
-// x0[k] -= deps; // differential step in coordinate k
-
-// Cell_Pass(1L, 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] =
-// (float) (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, double *x, vector *fjac, double *fvect)
-// {
-// int i, k;
-// long lastpos = 0L;
-// double x0[6], fx1[6], fx2[6];
-// const double deps = 1e-8; //stepsize for numerical differentiation
-
-// CopyVec(6L, x, x0);
-
-// Cell_Pass(1L, globval.Cell_nLoc, x0, &lastpos);
-// CopyVec((long) 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(1L, 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(1L, 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,float x[],int n,float tolx, float 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
-
-// ****************************************************************************/
-
-// #define FREERETURN {free_matrix(alpha,1,n,1,n);free_vector(bet,1,n); free_vector(fvect,1,n); free_ivector(indx,1,n);return;}
-
-// void Newton_RaphsonS(int ntrial, float x[], int n, float tolx)
-// {
-// int k, i, *indx, *ivector();
-// float errx, d, *bet, *fvect, **alpha, *vector(), **matrix();
-// void usrfun(), ludcmp(), lubksb(), free_ivector(), free_vector(),
-// free_matrix();
-
-// errx = (float) 0.0;
-// // NR arrays start from 1 and not 0 !!!
-// indx = ivector(1, n);
-// bet = vector(1, n);
-// fvect = vector(1, n);
-// alpha = matrix(1, n, 1, n);
-
-// for (k = 1; k <= ntrial; k++)
-// { // loop over number of iterations
-// computeFandJS(x, n, alpha, fvect); // supply function values at x in fvect and Jacobian matrix in fjac
-
-// // Jacobian -Id
-// for (i = 1; i <= n; i++)
-// alpha[i][i] -= (float) 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
-// ludcmp(alpha, n, indx, &d);
-// lubksb(alpha, n, indx, bet);
-// errx = (float) 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;
-// FREERETURN}
-// }
-// // check whever closed orbit found out
-// if ((k >= ntrial) && (errx >= tolx * 100))
-// {
-// status.codflag = false;
-// FREERETURN}
-// }
-
-// #undef FREERETURN
-
-// /****************************************************************************/
-// /* 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, double x[], int ntrial, double tolx)
-// {
-// int k, i;
-// double errx, bet[6], fvect[6];
-// vector alpha[6];
-
-// errx = 0.0;
-
-// for (k = 1; k <= ntrial; k++)
-// { // loop over number of iterations
-// computeFandJ(n,x, alpha, fvect); // supply function values at x in fvect and Jacobian matrix in fjac
-
-// // 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 Get_NAFF(int nterm, long ndata, double T[DIM][NTURN],
-// double *fx, double *fz, int nb_freq[2])
-
-// Purpose:
-// Compute quasiperiodic approximation of a phase space trajectory
-// using NAFF Algorithm ((c) Laskar, IMCCE)
-
-// Input:
-// nterm number of frequencies to look for
-// if not multiple of 6, truncated to lower value
-// ndata size of the data to analyse
-// T 6D vector to analyse
-
-// Output:
-// fx frequencies found in the H-plane
-// fz frequencies found in the V-plane
-// nb_freq number of frequencies found out in each plane
-
-// Return:
-// none
-
-// Global variables:
-// g_NAFVariable see modnaff.c
-// M_2_PI defined in math.h
-// trace ON or TRUE for debugging
-
-// Specific functions:
-// naf_initnaf, naf_cleannaf
-// naf_mftnaf
-
-// Comments:
-// none
-
-// ****************************************************************************/
-// /* Frequency Map Analysis */
-// #include "modnaff.h"
-// #include "complexeheader.h"
-// /* Analyse en Frequence */
-// double pi = M_PI;
-// void Get_NAFF(int nterm, long ndata, double T[DIM][NTURN],
-// double *fx, double *fz, int nb_freq[2])
-// {
-// /* Test whether ndata is divisible by 6 -- for NAFF -- */
-// /* Otherwise truncate ndata to lower value */
-// long r = 0; /* remainder of the euclidian division of ndata by 6 */
-// int i;
-
-// if ((r = ndata % 6) != 0) {
-// printf("Get_NAFF: Warning ndata = %ld, \n", ndata);
-// ndata -= r;
-// printf("New value for NAFF ndata = %ld \n", ndata);
-// }
-
-// g_NAFVariable.DTOUR = M_2_PI; /* size of one "cadran" */
-// g_NAFVariable.XH = M_2_PI; /* step */
-// g_NAFVariable.T0 = 0.0; /* time t0 */
-// g_NAFVariable.NTERM = nterm; /* max term to find */
-// g_NAFVariable.KTABS = ndata; /* number of data: must be a multiple of 6 */
-// g_NAFVariable.m_pListFen = NULL; /* no window */
-// g_NAFVariable.TFS = NULL; /* will contain frequency */
-// g_NAFVariable.ZAMP = NULL; /* will contain amplitude */
-// g_NAFVariable.ZTABS = NULL; /* will contain data to analyze */
-
-// /****************************************************/
-// /* internal use in naf */
-// g_NAFVariable.NERROR = 0;
-// g_NAFVariable.ICPLX = 1;
-// g_NAFVariable.IPRT = -1; /* 1 for diagnostics */
-// g_NAFVariable.NFPRT = stdout; /* NULL */
-// g_NAFVariable.NFS = 0;
-// g_NAFVariable.IW = 1;
-// g_NAFVariable.ISEC = 1;
-// g_NAFVariable.EPSM = 0;
-// g_NAFVariable.UNIANG = 0;
-// g_NAFVariable.FREFON = 0;
-// g_NAFVariable.ZALP = NULL;
-// g_NAFVariable.m_iNbLineToIgnore = 1; /* unused */
-// g_NAFVariable.m_dneps = 1.e10;
-// g_NAFVariable.m_bFSTAB = FALSE; /* unused */
-// /* end of interl use in naf */
-// /****************************************************/
-
-// /* NAFF initialization */
-// naf_initnaf();
-
-// /**********************/
-// /* Analyse in H-plane */
-// /**********************/
-
-// /* fills up complexe vector for NAFF analysis */
-// for(i = 0; i < ndata; i++) {
-// g_NAFVariable.ZTABS[i].reel = T[0][i]; /* x */
-// g_NAFVariable.ZTABS[i].imag = T[1][i]; /* xp */
-// }
-
-// /* Get out the mean value */
-// naf_smoy(g_NAFVariable.ZTABS);
-// naf_prtabs(g_NAFVariable.KTABS,g_NAFVariable.ZTABS, 20);
-// // trace=on;
-// naf_mftnaf(nterm,fabs(g_NAFVariable.FREFON)/g_NAFVariable.m_dneps);
-
-// /* fill up H-frequency vector */
-// for (i = 1; i <= g_NAFVariable.NFS; i++) {
-// fx[i-1] = g_NAFVariable.TFS[i];
-// }
-
-// nb_freq[0] = g_NAFVariable.NFS; /* nb of frequencies found out by NAFF */
-
-// if (trace) /* print out results */
-// {
-// printf("(x,x') phase space: NFS=%d\n",g_NAFVariable.NFS);
-// for (i = 1; i <= g_NAFVariable.NFS; i++) {
-// printf("AMPL=%15.8E+i*%15.8E abs(AMPL)=%15.8E arg(AMPL)=%15.8E FREQ=%15.8E\n",
-// g_NAFVariable.ZAMP[i].reel,g_NAFVariable.ZAMP[i].imag,
-// i_compl_module(g_NAFVariable.ZAMP[i]),
-// i_compl_angle(g_NAFVariable.ZAMP[i]),
-// g_NAFVariable.TFS[i]);
-// }
-// }
-
-// /**********************/
-// /* Analyse in V-plane */
-// /**********************/
-
-// /* fill up complexe vector for NAFF analysis */
-// for (i = 0; i < ndata; i++) {
-// g_NAFVariable.ZTABS[i].reel = T[2][i]; /* z */
-// g_NAFVariable.ZTABS[i].imag = T[3][i]; /*zp */
-// }
-
-// naf_mftnaf(nterm,fabs(g_NAFVariable.FREFON)/g_NAFVariable.m_dneps);
-
-// /* fills up V-frequency vector */
-// for (i = 1; i <= g_NAFVariable.NFS; i++) {
-// fz[i-1] = g_NAFVariable.TFS[i];
-// }
-
-// nb_freq[1] = g_NAFVariable.NFS; /* nb of frequencies found out by NAFF */
-
-// if (trace) /* print out results */
-// {
-// printf("(z,z') phase space: NFS=%d\n",g_NAFVariable.NFS);
-// for (i = 1; i <= g_NAFVariable.NFS; i++) {
-// printf("AMPL=%15.8E+i*%15.8E abs(AMPL)=%15.8E arg(AMPL)=%15.8E FREQ=%15.8E\n",
-// g_NAFVariable.ZAMP[i].reel,g_NAFVariable.ZAMP[i].imag,
-// i_compl_module(g_NAFVariable.ZAMP[i]),
-// i_compl_angle(g_NAFVariable.ZAMP[i]),
-// g_NAFVariable.TFS[i]);
-// }
-// }
-
-// /* free out memory used by NAFF */
-// naf_cleannaf();
-// }
-
-// /****************************************************************************/
-// /* void Get_Tabshift(double Tab[DIM][NTURN], double Tab0[DIM][NTURN],
-// long nbturn, long nshift)
-
-// Purpose: used by fmap
-// Store in Tab0 values of Tab shifted by nshift turn
-
-// Input:
-// Tab tracking tabular
-// nshift nb of turns to shift
-// nbturn nb of turns
-
-// Output:
-// Tab0 tracking tabular
-
-// Return:
-// none
-
-// Global variables:
-// none
-
-// Specific functions:
-// none
-
-// Comments:
-
-
-// ****************************************************************************/
-// void Get_Tabshift(double Tab[DIM][NTURN], double Tab0[DIM][NTURN], long nbturn, long nshift)
-// {
-// long k = 0L, k1 = 0L;
-
-// for (k = 0L; k < nbturn ; k++){
-// k1 = k + nshift;
-// Tab0[0][k] = Tab[0][k1];
-// Tab0[1][k] = Tab[1][k1];
-// Tab0[2][k] = Tab[2][k1];
-// Tab0[3][k] = Tab[3][k1];
-// }
-
-// }
-
-// /****************************************************************************/
-// /* void Get_freq(double *fx, double *fz, double *nux, double *nuz)
-
-// Purpose: used by fmap
-// Looks for tunes from frequency vectors output from NAFF
-
-// Input:
-// fx vector of horizontal frequencies
-// fz vector of verticcal frequencies
-
-// Output:
-// nux H tune
-// nuz V tune
-
-// Return:
-// none
-
-// Global variables:
-// none
-
-// Specific functions:
-// none
-
-// Comments:
-
-
-// ****************************************************************************/
-// void Get_freq(double *fx, double *fz, double *nux, double *nuz)
-// {
-// const double eps0 = 1e-4;
-// const double eps1 = 1e-6;
-
-// // case of nux
-// if (fabs(fx[0]) < eps0){
-// *nux = fx[1];
-// }
-// else {
-// *nux = fx[0];
-// }
-
-// // case of nuz
-// if (fabs(fz[0]) < eps0) {
-// if (fabs(fabs(fz[1]) - fabs(*nux)) < eps1) {
-// if (fabs(fabs(fz[2]) - fabs(*nux)) < eps1) {
-// *nuz = fz[3];
-// }
-// else *nuz = fz[2];
-// }
-// else *nuz = fz[1];
-// }
-// else{
-// if (fabs(fabs(fz[0]) - fabs(*nux)) < eps1) {
-// if (fabs(fabs(fz[1]) - fabs(*nux)) < eps1) {
-// *nuz = fz[2];
-// }
-// else *nuz = fz[1];
-// }
-// else *nuz = fz[0];
-// }
-// *nuz = fabs(*nuz); *nux = fabs(*nux);
-// }
-
-
-// /* DO NOT REMOVE FOR f2c ...*/
-// /* To be fixed later .. */
-// void MAIN__(void)
-// {
-
-// }
/****************************************************************************
*
diff --git a/tracy/tracy/src/soleillib.cc b/tracy/tracy/src/soleillib.cc
index a85cab4..fe1d051 100644
--- a/tracy/tracy/src/soleillib.cc
+++ b/tracy/tracy/src/soleillib.cc
@@ -42,6 +42,7 @@ void Get_Disp_dp(void)
FILE *outf;
double dP = 0e0;
CellType Cell;
+ long lastpos =0L;
if (trace) fprintf(stdout,"Entering Get_Disp_dp function ...\n");
@@ -53,8 +54,7 @@ void Get_Disp_dp(void)
for (i = 1; i <= 20; i++) {
dP = -0.003 + 1e-6 + i*0.0006;
- //~ getcod(dP, lastpos);
- findcod(dP);
+ getcod(dP, lastpos);
Ring_GetTwiss(true, dP); /* Compute and get Twiss parameters */
getelem(0, &Cell);
fprintf(outf,"%+e %+e %+e\n", dP, Cell.BeamPos[0], Cell.Eta[0]);
@@ -251,10 +251,10 @@ void Hcofonction(long pos, double dP)
{
CellType Cell;
long i;
- long lastpos = 1;
+ long lastpos = 1L;
Vector2 H;
- //~ getcod(dP, lastpos); /* determine closed orbit */
- findcod(dP);
+ getcod(dP, lastpos); /* determine closed orbit */
+
if (lastpos != globval.Cell_nLoc) printf("Ring unstable for dp=%+e @ pos=%ld\n", dP, lastpos);
Ring_GetTwiss(true, dP); /* Compute and get Twiss parameters */
@@ -1328,8 +1328,8 @@ void NuDp(long Nb, long Nbtour, double emax)
status = true;
}
- //~ getcod(dp, lastpos); // get cod for printout
- findcod(dp);
+ long lastpos=0L;
+ getcod(dp, lastpos); // get cod for printout
fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
@@ -1709,8 +1709,7 @@ void Enveloppe(double x, double px, double y, double py, double dp, double nturn
CellType Cell;
/* Get cod the delta = energy*/
- //~ getcod(dp, lastpos);
- findcod(dp);
+ getcod(dp, lastpos);
printf("xcod=%.5e mm zcod=% .5e mm \n", globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
@@ -3266,8 +3265,8 @@ void MomentumAcceptance(long deb, long fin, double ep_min, double ep_max,
do
{
- //~ getcod(dP=0.0, lastpos); /* determine closed orbit */
- findcod(dP=0.0);
+ getcod(dP=0.0, lastpos); /* determine closed orbit */
+
getelem(pos,&Cell);
// coordinates around closed orbit which is non zero for 6D tracking
x = Cell.BeamPos[0];
@@ -3419,8 +3418,8 @@ void MomentumAcceptance(long deb, long fin, double ep_min, double ep_max,
/***************************************************************/
do {
- //~ getcod(dP=0.0, lastpos); /* determine closed orbit */
- findcod(dP=0.0);
+ getcod(dP=0.0, lastpos); /* determine closed orbit */
+
getelem(pos,&Cell);
// coordinates around closed orbit which is non zero for 6D tracking
x = Cell.BeamPos[0];
@@ -3505,14 +3504,14 @@ void MomentumAcceptance(long deb, long fin, double ep_min, double ep_max,
****************************************************************************/
void set_vectorcod(Vector codvector[], double dP)
{
- long k = 0;
+ long k = 0L, lastpos = 0L;
CellType Cell;
Vector zerovector;
zerovector.zero();
- //~ getcod(dP, lastpos); /* determine closed orbit */
- findcod(dP);
+ getcod(dP, lastpos); /* determine closed orbit */
+
if (status.codflag == 1) { /* cod exists */
for (k = 1L; k <= globval.Cell_nLoc; k++){
--
GitLab