From 7009f839e896d80ed93e34277688bf79b6c71192 Mon Sep 17 00:00:00 2001
From: nadolski <nadolski@9a6e40ed-f3a0-4838-9b4a-bf418f78e88d>
Date: Fri, 1 Oct 2010 16:44:49 +0000
Subject: [PATCH] Ring_GetChroma modify. As in Tracy II, convergence test
 introduced since at commonDP (10-10) the computing precision is low.

---
 tracy/tracy/src/t2ring.cc | 137 +++++++++++++++++++++++++-------------
 1 file changed, 92 insertions(+), 45 deletions(-)

diff --git a/tracy/tracy/src/t2ring.cc b/tracy/tracy/src/t2ring.cc
index 00b0a46..b88a96d 100644
--- a/tracy/tracy/src/t2ring.cc
+++ b/tracy/tracy/src/t2ring.cc
@@ -475,56 +475,103 @@ void Cell_Twiss(long i0, long i1, ss_vect<tps> &Ascr, bool chroma, bool ring,
        WARNING : this routine does not give the chromaticities for dP != 0
                    but the local slope of the curve xi=f(dP)                   
        16/10/03 Modified convergence test: now done for both planes
+       01/09/10 Modify the convergence criteria on relative diff of the chroma value
+                The previous test does not work well for non zero chromaticities
 ****************************************************************************/
 #define n               4
-void Ring_Getchrom(double dP)
-{
-  long int  lastpos = 0;
-  int       j;
-  Vector2   alpha = {0.0, 0.0}, beta = {0.0, 0.0}, gamma = {0.0, 0.0};
-  Vector2   nu = {0.0, 0.0}, nu0 = {0.0, 0.0};
-
-  if (dP != 0.0)
-    fprintf(stdout,"Ring_Getchrom: Warning this is NOT the CHROMA, dP=%e\n",
-	    dP);
-  
-  /* Get cod for energy dP - globval.dPcommon*/
-  GetCOD(globval.CODimax, globval.CODeps, dP-globval.dPcommon*0.5, lastpos);
-  
-  if (!status.codflag) {
-    /* if no cod */
-    fprintf(stdout,"Ring_Getchrom:  Lattice is unstable for"
-	    " dP-globval.dPcommon=% .5e\n", dP-globval.dPcommon*0.5);
-    return;
-  }
-  
-  /* get tunes for energy dP - globval.dPcommon/2 from oneturn matrix */
-  Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu0);
-  
-  /* Get cod for energy dP+globval.dPcommon*/
-  GetCOD(globval.CODimax, globval.CODeps, dP+globval.dPcommon*0.5, lastpos);
-  
-  if (!status.codflag) { /* if no cod */
-    fprintf(stdout,"Ring_Getchrom  Lattice is unstable for"
-	    " dP+globval.dPcommon=% .5e \n", dP+globval.dPcommon*0.5);
-    return;
-  }
-
-  /* get tunes for energy dP+globval.dPcommon/2 from one turn matrix */
-  Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu);
-  
-  if (!globval.stable) {
-    printf("Ring_Getchrom:  Lattice is unstable\n");
-  }
-
-  /* Get chromaticities by numerical differentiation*/
-  for (j = 0; j <= 1; j++)
-    globval.Chrom[j] = (nu[j]-nu0[j])/globval.dPcommon;
-  
-  status.chromflag = true;
+#define chromeps        1e-6 /* convergence condition for chromaticity computation */
+#define LOG10 log(10.0)
+
+void Ring_Getchrom(double dP) {
+	long int lastpos = 0;
+	double NORMchro = 1.0;
+	double dPlocal = 0.0, expo = 0.0, Norm_TEMP = 0.0;
+	int j;
+	Vector2 alpha = { 0.0, 0.0 }, beta = { 0.0, 0.0 }, gamma = { 0.0, 0.0 };
+	Vector2 nu = { 0.0, 0.0 }, nu0 = { 0.0, 0.0 }, TEMP={0.0,0.0}, Chrom={0.0,0.0};;
+
+	if (dP != 0.0)
+		fprintf(stdout,
+				"Ring_Getchrom: Warning this is NOT the CHROMA, dP=%e\n", dP);
+
+	/* initialization */
+	globval.Chrom[0] = 1e38;
+	globval.Chrom[1] = 1e38;
+
+	expo = log(globval.dPcommon) / LOG10;
+	do {
+		for (j = 0; j <= 1; j++)
+			Chrom[j] = globval.Chrom[j];
+
+		dPlocal = exp(expo * LOG10);
+		//Get cod for energy dP - globval.dPcommon
+		GetCOD(globval.CODimax, globval.CODeps, dP - dPlocal * 0.5, lastpos);
+
+		if (!status.codflag) {
+			// if no cod
+			fprintf(stdout, "Ring_Getchrom:  Lattice is unstable for"
+				" dP-globval.dPcommon=% .5e\n", dP - dPlocal * 0.5);
+			return;
+		}
+
+		//get tunes for energy dP - globval.dPcommon/2 from oneturn matrix
+		Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu0);
+
+		//Get cod for energy dP+globval.dPcommon
+		GetCOD(globval.CODimax, globval.CODeps, dP + dPlocal * 0.5, lastpos);
+
+		if (!status.codflag) { // if no cod
+			fprintf(stdout, "Ring_Getchrom  Lattice is unstable for"
+				" dP+globval.dPcommon=% .5e \n", dP + dPlocal * 0.5);
+			return;
+		}
+
+		//get tunes for energy dP+globval.dPcommon/2 from one turn matrix
+		Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu);
+
+		if (!globval.stable) {
+			printf("Ring_Getchrom:  Lattice is unstable\n");
+		}
+
+		//Get chromaticities by numerical differentiation
+		for (j = 0; j <= 1; j++)
+			globval.Chrom[j] = (nu[j] - nu0[j]) / globval.dPcommon;
+
+		/* Get chromaticities by numerical differentiation*/
+		for (j = 0; j <= 1; j++) {
+			globval.Chrom[j] = (nu[j] - nu0[j]) / dPlocal;
+			TEMP[j] = fabs(globval.Chrom[j] - Chrom[j]);
+		}
+
+		Norm_TEMP = sqrt(TEMP[0] * TEMP[0] + TEMP[1] * TEMP[1]);
+		NORMchro = sqrt(globval.Chrom[0] * globval.Chrom[0] + globval.Chrom[1]
+				* globval.Chrom[1]);
+
+		// TEST CHROMA convergence
+		if (trace) {
+			fprintf(
+					stdout,
+					"\nexpo % e xix = % e xiz = % e Norm_TEMP = %e TEMPX %+e TEMPZ %+e\n",
+					expo, Chrom[0], Chrom[1], Norm_TEMP, TEMP[0], TEMP[1]);
+			fprintf(stdout, "expo % e nux = % e nuz = % e dPlocal= %+e\n",
+					expo, nu0[0], nu0[1], dP - 0.5 * dPlocal);
+			fprintf(stdout, "expo % e nux = % e nuz = % e dPlocal= %+e\n",
+					expo, nu[0], nu[1], dP + 0.5 * dPlocal);
+		}
+		expo += 1.0;
+		if (trace)
+			fprintf(stdout, "%+e %+.12e %+.12e %+.12e %+.12e\n", dPlocal,
+					globval.Chrom[0], fabs(globval.Chrom[0] - Chrom[0])
+							/ Chrom[0], globval.Chrom[1], fabs(globval.Chrom[1]
+							- Chrom[1]) / Chrom[1]);
+	} while ((Norm_TEMP > chromeps * NORMchro) || (Norm_TEMP == 0e0));
+
+	status.chromflag = true;
 }
 
 #undef n
+#undef chromeps
+#undef LOG10
 
 
 /****************************************************************************/
-- 
GitLab