From 29b7e973e2baaa29d5512af6675ebf3bf50ff550 Mon Sep 17 00:00:00 2001
From: nadolski <nadolski@9a6e40ed-f3a0-4838-9b4a-bf418f78e88d>
Date: Tue, 10 Feb 2015 15:18:50 +0000
Subject: [PATCH] introduction of asym aperture validated

---
 tracy/tracy/inc/tracy_global.h |   1 +
 tracy/tracy/src/prtmfile.cc    |   5 +-
 tracy/tracy/src/rdmfile.cc     |   2 +
 tracy/tracy/src/t2cell.cc      |  13 +++
 tracy/tracy/src/t2elem.cc      | 171 +++++++++++++++++++++++++++------
 tracy/tracy/src/t2lat.cc       | 111 ++++++++++++++++++---
 6 files changed, 258 insertions(+), 45 deletions(-)

diff --git a/tracy/tracy/inc/tracy_global.h b/tracy/tracy/inc/tracy_global.h
index 8e85819..812b133 100644
--- a/tracy/tracy/inc/tracy_global.h
+++ b/tracy/tracy/inc/tracy_global.h
@@ -76,6 +76,7 @@ struct DriftType {
 
 struct ApertureType {
   double Aper[4][3]; // Aperture
+  double LSN;
 };
 
 
diff --git a/tracy/tracy/src/prtmfile.cc b/tracy/tracy/src/prtmfile.cc
index d13439b..96b10ec 100644
--- a/tracy/tracy/src/prtmfile.cc
+++ b/tracy/tracy/src/prtmfile.cc
@@ -87,7 +87,7 @@
  type             predefined integer number for different type of element
  Type codes:
  marker         -1
- drift	    0
+ drift	         0
  multipole       1
  cavity          2
  thin kick       3
@@ -171,7 +171,7 @@ void prtAper(FILE *fp, const double (*Aper)[3]) {
   int i = 0;
 
   for (i = 0; i <= 3; i++) {
-      fprintf(fp, "%23.16e %23.16e %23.16e\n", Aper[i][0], Aper[i][1],
+    fprintf(fp, "%+23.16e %+23.16e %+23.16e\n", Aper[i][0], Aper[i][1],
       Aper[i][2]);
   }
 }
@@ -246,6 +246,7 @@ void prtmfile(const char mfile_dat[]) {
       break;
     case Aperture:
       prtName(mfile, i, aperture_, 0, 0);
+      // fprintf(mfile, " L: %23.16e\n", Cell[i].Elem.PL);
       prtAper(mfile, Cell[i].Elem.A->Aper);
       break;
     case Insertion:
diff --git a/tracy/tracy/src/rdmfile.cc b/tracy/tracy/src/rdmfile.cc
index 11c52d4..9d60527 100644
--- a/tracy/tracy/src/rdmfile.cc
+++ b/tracy/tracy/src/rdmfile.cc
@@ -187,6 +187,8 @@ void rdmfile(const char *mfile_dat)
             break;
         case marker:
             break;
+        case Aperture:
+            break;
         case drift:
             inf.getline(line, max_str);
             if (prt)
diff --git a/tracy/tracy/src/t2cell.cc b/tracy/tracy/src/t2cell.cc
index b334d74..4989750 100644
--- a/tracy/tracy/src/t2cell.cc
+++ b/tracy/tracy/src/t2cell.cc
@@ -92,6 +92,9 @@ void Elem_Pass(const long i, ss_vect<T> &x)
     case marker:
       Marker_Pass(Cell[i], x);
       break;
+    case Aperture:
+      Aperture_Pass(Cell[i], x);
+      break;
     case Spreader:
       break;
     case Recombiner:
@@ -136,6 +139,8 @@ void Elem_Pass_M(const long i, Vector &xref, Matrix &x)
       break;
     case marker:   /* nothing */
       break;
+    case Aperture:   /* nothing */
+      break;
     default:
       fprintf(stdout,"Elem_Pass_M: ** undefined type\n");
       break;
@@ -178,6 +183,8 @@ void Cell_SetdP(const double dP)
       break;
     case marker:
       break;
+    case Aperture:
+      break;
     case Spreader:
       break;
     case Recombiner:
@@ -411,6 +418,9 @@ bool linearel(long i)
   case Cavity:
     status = true;
     break;
+  case Aperture:
+    status = true;
+    break;
   case marker:
     status = true;
     break;
@@ -536,6 +546,9 @@ void Cell_Concat(double dP)
         case marker:   /* nothing */
           break;
 
+        case Aperture:   /* nothing */
+          break;
+
         default:
           printf("**Cell_Concat: undefined type\n");
 	  break;
diff --git a/tracy/tracy/src/t2elem.cc b/tracy/tracy/src/t2elem.cc
index 6de56a2..2cfddda 100644
--- a/tracy/tracy/src/t2elem.cc
+++ b/tracy/tracy/src/t2elem.cc
@@ -1187,36 +1187,65 @@ void Aperture_Pass(CellType &Cell, ss_vect<T> &X) {
     ApertureType *A;
     bool outofrange = false;
 
-    bool trace = true;
+    bool trace = false;
     
     elemp = &Cell.Elem;
     A = elemp->A;
 
     /* Global -> Local */
     GtoL(X, Cell.dS, Cell.dT, 0.0, 0.0, 0.0);
+        
     
     /* horizontal plane aperture */
-    if (X[y_] < A->Aper[0][1] & X[y_] > A->Aper[0][2])
-    	if (X[x_] > A->Aper[0][0]) {outofrange = true;}
-    if (X[y_] < A->Aper[1][1] & X[y_] > A->Aper[1][2])
-    	if (X[x_] > A->Aper[1][0]) {outofrange = true;}
+    /* part lost if x>xplim and (z>zpp ou z< zpm) */
+    if (X[x_] > A->Aper[0][0]) 
+      if (X[y_] < A->Aper[0][1] || X[y_] > A->Aper[0][2])
+    	 {outofrange = true;}
+
+    /* part lost if x<xmlim and (z>zmp ou z< zmm) */
+    if (X[x_] < A->Aper[1][0]) 
+      if (X[y_] < A->Aper[1][1] || X[y_] > A->Aper[1][2])
+    	{outofrange = true;}
     
-    if (trace)
-    cout << "X[y_] " << X[y_] << endl;
+    if (outofrange && trace){ /* for debugging */
+    	fprintf(stdout,"\n x[m]=%+23.16e z[m]=%+23.16e\n", X[x_], X[y_]); 
+        fprintf(stdout, "OUT OF RANGE in V-plane \n");
+		fprintf(stdout,"xplim=%+23.16e zmp=%+23.16e zpp=%+23.16e\n", 
+		 A->Aper[0][0],  A->Aper[0][1],  A->Aper[0][2]) ;
+
+		fprintf(stdout,"xmlim=%+23.16e zmm=%+23.16e zmp=%+23.16e\n", 
+		 A->Aper[1][0],  A->Aper[1][1],  A->Aper[1][2]) ;
+	 }
     
     if (outofrange) {
-      X[x_] = NAN;
+      X[y_] = NAN;
       return;
     }
     
     /* vertical plane aperture */
-    if (X[x_] < A->Aper[2][1] & X[x_] > A->Aper[2][2])
-    	if (X[y_] > A->Aper[2][0]) {outofrange = true;}
-    if (X[x_] < A->Aper[3][1] & X[x_] > A->Aper[3][2])
-    	if (X[y_] > A->Aper[3][0]) {outofrange = true;}
+    /* part lost if z>zplim and (x>xpp ou x< xpm) */
+    if (X[y_] > A->Aper[2][0])  
+      if (X[x_] < A->Aper[2][1] || X[x_] > A->Aper[2][2])
+        {outofrange = true;}
+    /* part lost if z<zmlim and (x>xmp ou x< xmm) */
+    if (X[y_] < A->Aper[3][0])
+      if (X[x_] < A->Aper[3][1] || X[x_] > A->Aper[3][2])
+    	{outofrange = true;}
+ 
     
+    if (outofrange && trace){ /* for debugging */
+    	fprintf(stdout,"\n x[m]=%+23.16e z[m]=%+23.16e\n", X[x_], X[y_]); 
+        fprintf(stdout, "OUT OF RANGE in H-plane\n");
+
+		fprintf(stdout,"yplim=%+23.16e xpm=%+23.16e xpp=%+23.16e\n", 
+		 A->Aper[2][0],  A->Aper[2][1],  A->Aper[2][2]) ;
+
+		fprintf(stdout,"ymlim=%+23.16e xmm=%+23.16e xmp=%+23.16e\n", 
+		 A->Aper[3][0],  A->Aper[3][1],  A->Aper[1][2]) ;
+    }
+
     if (outofrange) {
-      X[y_] = NAN;
+      X[x_] = NAN;
       return;
     }
     
@@ -3222,6 +3251,34 @@ void SI_init(void) {
             * m_e)) * pow(globval.Energy, 5.0);
 }
 
+/****************************************************************************
+ * void Mpole_Print(FILE **f, long Fnum1)
+
+ Purpose: called by Elem_Print
+ Print out drift characteristics in a file
+ Name, kind, length, method, slice number
+
+ Input:
+ Fnum1 Family number
+ f     pointer on file id
+
+ Output:
+ none
+
+ Return:
+ none
+
+ Global variables:
+ none
+
+ Specific functions:
+ none
+
+ Comments:
+ none
+
+ ****************************************************************************/
+
 static void Mpole_Print(FILE *f, int Fnum1) {
     elemtype *elemp;
     MpoleType *M;
@@ -3234,6 +3291,51 @@ static void Mpole_Print(FILE *f, int Fnum1) {
     fprintf(f, "   Method: %d, N=%4d\n", M->Pmethod, M->PN);
 }
 
+
+/****************************************************************************
+ * void Mpole_Print(FILE **f, long Fnum1)
+
+ Purpose: called by Elem_Print
+ Print out drift characteristics in a file
+ Name, kind, length, method, slice number
+
+ Input:
+ Fnum1 Family number
+ f     pointer on file id
+
+ Output:
+ none
+
+ Return:
+ none
+
+ Global variables:
+ none
+
+ Specific functions:
+ none
+
+ Comments:
+ none
+
+ ****************************************************************************/
+static void Aperture_Print(FILE *f, int Fnum1) {
+    elemtype *elemp;
+    ApertureType *A;
+    long i = 0;
+
+    elemp = &ElemFam[Fnum1 - 1].ElemF;
+    A = elemp->A;
+    fprintf(f, "Element[%3d ] \n", Fnum1);
+    fprintf(f, "   Name: %.*s,  Kind:   mpole,  L=% .8E\n", SymbolLength,
+            elemp->PName, elemp->PL);
+    for (i = 0; i <= 3; i++) {
+      fprintf(f, "%+23.16e %+23.16e %+23.16e\n", A->Aper[i][0], A->Aper[i][1],
+        A->Aper[i][2]);
+    }
+}
+
+
 /****************************************************************************
  * void Drift_Print(FILE **f, long Fnum1)
 
@@ -3391,13 +3493,15 @@ void Elem_Print(FILE *f, int Fnum1) {
     case drift:
         Drift_Print(f, Fnum1);
         break;
-
     case Mpole:
         Mpole_Print(f, Fnum1);
         break;
     case Wigl:
         Wiggler_Print(f, Fnum1);
         break;
+    case Aperture:
+        Aperture_Print(f, Fnum1);
+        break;
     case FieldMap:
         break;
     case Insertion:
@@ -3463,6 +3567,9 @@ double Elem_GetKval(int Fnum1, int Knum1, int Order) {
         case marker:
             Result = 0.0;
             break;
+        case Aperture:
+            Result = 0.0;
+            break;
         case Cavity:
             Result = 0.0;
             break;
@@ -3591,27 +3698,28 @@ void Drift_Alloc(elemtype *Elem) {
  ****************************************************************************/
  void Aperture_Alloc(elemtype *Elem) {
 
-    ApertureType *A;;
+    ApertureType *WITH8;
 
     /* Memory allocation */
     Elem->A = (ApertureType *) malloc(sizeof(ApertureType));
-    A = Elem->A;
+    WITH8 = Elem->A;
     
-    A->Aper[0][0] = +1.0;
-    A->Aper[0][1] = -1.0;
-    A->Aper[0][2] = +1.0;
-
-    A->Aper[1][0] = -1.0;
-    A->Aper[1][1] = -1.0;
-    A->Aper[1][2] = +1.0;
-
-    A->Aper[2][0] = +1.0;
-    A->Aper[2][1] = -1.0;
-    A->Aper[2][2] = +1.0;
+    WITH8->LSN = -1.0;
+    WITH8->Aper[0][0] = +1.0;
+    WITH8->Aper[0][1] = -1.0;
+    WITH8->Aper[0][2] = +1.0;
+
+    WITH8->Aper[1][0] = -1.0;
+    WITH8->Aper[1][1] = -1.0;
+    WITH8->Aper[1][2] = +1.0;
+
+    WITH8->Aper[2][0] = +1.0;
+    WITH8->Aper[2][1] = -1.0;
+    WITH8->Aper[2][2] = +1.0;
       
-    A->Aper[3][0] = -1.0;
-    A->Aper[3][1] = -1.0;
-    A->Aper[3][2] = +1.0;
+    WITH8->Aper[3][0] = -1.0;
+    WITH8->Aper[3][1] = -1.0;
+    WITH8->Aper[3][2] = +1.0;
 
 }
 
@@ -4302,7 +4410,8 @@ void Aperture_Init(int Fnum1) {
         /* Memory allocation and set everything to zero */
         Aperture_Alloc(&cellp->Elem);
         memcpy(cellp->Elem.PName, elemfamp->ElemF.PName, sizeof(partsName));
-
+        *cellp->Elem.A = *elemfamp->ElemF.A;
+        
         cellp->dT[0] = 1.0;
         cellp->dT[1] = 0.0;
         cellp->dS[0] = 0.0;
diff --git a/tracy/tracy/src/t2lat.cc b/tracy/tracy/src/t2lat.cc
index 4c76fd2..5a11da3 100644
--- a/tracy/tracy/src/t2lat.cc
+++ b/tracy/tracy/src/t2lat.cc
@@ -1715,6 +1715,48 @@ static bool Lat_CheckWiggler(FILE **fo, long i, struct LOC_Lattice_Read *LINK)
   return true;
 }
 
+static bool Lat_CheckAperture(FILE **fo, long i, struct LOC_Lattice_Read *LINK)
+{
+  bool         Result;
+  int j;
+  ElemFamType  *WITH;
+  elemtype     *WITH1;
+  ApertureType  *WITH8;
+
+  Result = false;
+  WITH = &ElemFam[i-1]; WITH1 = &WITH->ElemF; WITH8 = WITH1->A;
+   
+  if (trace) 
+  for (j = 0; j <= 3; j++)  {
+    printf("Lat_CheckAperture: %+f %f %+f \n", WITH8->Aper[j][0], 
+                              WITH8->Aper[j][1], WITH8->Aper[j][2]);
+  }
+
+  if (WITH8->Aper[0][0]>= WITH8->Aper[1][0] &&
+      WITH8->Aper[2][0]>= WITH8->Aper[3][0] &&
+      WITH8->Aper[0][1]<= WITH8->Aper[0][2] &&  
+      WITH8->Aper[1][1]<= WITH8->Aper[1][2] && 
+      WITH8->Aper[2][1]<= WITH8->Aper[2][2]  && 
+      WITH8->Aper[3][1]<= WITH8->Aper[3][2]) return true;
+  printf("\n");
+  printf(">>> Incorrect definition of %.*s\n\n", NameLength, WITH1->PName);
+  printf("    Asymmetric absorber : check limits \n");
+  if (WITH8->Aper[0][0]< WITH8->Aper[1][0])
+      printf("warning Aper[0][0] < Aper[1][0] \n");
+  if (WITH8->Aper[2][0]< WITH8->Aper[3][0])
+      printf("warning Aper[2][0] < Aper[3][0] \n");
+  if (WITH8->Aper[0][1]> WITH8->Aper[0][2])
+      printf("warning Aper[0][1] > Aper[0][2] \n");
+  if (WITH8->Aper[1][1]> WITH8->Aper[1][2])
+      printf("warning Aper[1][1] > Aper[1][2] \n");
+  if (WITH8->Aper[2][1]> WITH8->Aper[2][2])
+      printf("warning Aper[2][1] > Aper[2][2] \n");
+  if (WITH8->Aper[3][1]> WITH8->Aper[3][2])
+      printf("warning Aper[3][1] > Aper[3][2] \n");
+  exit_(1);
+  return true;
+}
+
 /* Local variables for Lat_DealElement: */
 struct LOC_Lat_DealElement
 {
@@ -1801,6 +1843,12 @@ static void CheckWiggler(long i, struct LOC_Lat_DealElement *LINK)
     longjmp(LINK->_JL9999, 1);
 }
 
+static void CheckAperture(long i, struct LOC_Lat_DealElement *LINK)
+{
+  if (!Lat_CheckAperture(LINK->fo, i, LINK->LINK))
+    longjmp(LINK->_JL9999, 1);
+}
+
 static void GetHOM(struct LOC_Lat_DealElement *LINK)
 {
   long    i;
@@ -1828,8 +1876,7 @@ static void GetHOM(struct LOC_Lat_DealElement *LINK)
 static void GetAPER(struct LOC_Lat_DealElement *LINK)
 {
   long i = 0L;
-  double  lim_;
-  double  pos_, neg_;
+  double  lim_ = 0.0, pos_ = 0.0, neg_=0.0;
   symset  SET;
 
   getest__(P_expset(SET, 1 << ((long)lparent)), "<(> expected", LINK);
@@ -1855,7 +1902,7 @@ static void ClearHOMandDBN(struct LOC_Lat_DealElement *LINK)
   long i;
 
   for (i = -HOMmax; i <= HOMmax; i++) {
-    LINK->B[i + HOMmax] = 0.0;
+    LINK->B[i + HOMmax]  = 0.0;
     LINK->BA[i + HOMmax] = false;
   }
   LINK->DBNsavemax = 0;
@@ -1876,6 +1923,8 @@ static void AssignHOM(long elem, struct LOC_Lat_DealElement *LINK)
   }
 }
 
+/* Set asymmetric aperture read from lattice file */
+
 static void AssignAPER(long elem, struct LOC_Lat_DealElement *LINK)
 {
   long       i;
@@ -1886,7 +1935,8 @@ static void AssignAPER(long elem, struct LOC_Lat_DealElement *LINK)
     A->Aper[i][0] = LINK->Aper[i][0]; 
     A->Aper[i][1] = LINK->Aper[i][1]; 
     A->Aper[i][2] = LINK->Aper[i][2]; 
-    if (!trace) printf("Aper[%ld] = %lf\n", i, A->Aper[i][0]);   
+    if (trace) printf("AssignAPER: %+f %f %+f \n", A->Aper[i][0], 
+                              A->Aper[i][1],A->Aper[i][2]);
   }
 }
 
@@ -2081,6 +2131,7 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
   InsertionType  *WITH5;   /* ID LSN */
   SolenoidType   *WITH7;
   ApertureType   *WITH8;   /* Asymmetric aperture LSN */
+  double Aper_[4][3]; // Aperture
   char str1[100] = "";
   char str2[100] = "";
   bool firstflag  = false; // flag for first kick input
@@ -2721,6 +2772,7 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
 
       test__(P_expset(SET, 1 << ((long)semicolon)), "<;> expected.", &V);
     }
+    
     GetSym__(&V);
     globval.Elem_nFam++;
     if (globval.Elem_nFam <= Elem_nFamMax) {
@@ -2837,7 +2889,7 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
 
     Example
 
-      ABS: aperture, aper=(+0.01, -0.01,  0.01,  zplim, xpm, xpp
+      ABS: aperture, aper=( +1.00, -1.00, +1.00,   xplim, zpm, zpp
                             -0.02, -0.01, -0.01,  zmlim, xmm, xmp
                             +0.02, -0.01, -0.01,  xplim, zpm, zpp
                             -0.02, -0.01, -0.01); xmlim, zmm, zmp
@@ -2846,27 +2898,62 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
     **************************************************************************/
 
   case aprsym:
+    
+    QL = 0e0; /* Length */
+ 
     getest__(P_expset(SET, 1 << ((long)comma)), "<comma> expected", &V);
-    getest__(P_addset(P_expset(SET1, 0), (long)apersym), "<aper> expected", &V);
-    getest__(P_expset(SET, 1 << ((long)eql)), "<=> expected", &V);
-    GetAPER(&V); /* read aperture */
+    
+    if (*V.sym == comma){
+        GetSym__(&V);
+        P_addset(P_expset(mysys, 0), (long)lsym);
+        P_addset(mysys, (long)apersym);
+		do {
+		  test__(mysys, "illegal parameter", &V); 
+		  sym1 = *V.sym;
+		  getest__(P_expset(SET, 1 << ((long)eql)), "<=> expected", &V);
+	  
+		  switch (sym1) {   /*11*/
+
+		  case lsym:
+			 QL = EVAL_(&V); /* get length */
+			 break;
+
+		  case apersym:
+			GetAPER(&V); /* read aperture */
+		  break;  
+	  
+		 default:
+		 break;
+		 }
+		 test__(P_expset(SET, (1 << ((long)comma)) | (1 << ((long)semicolon))),
+			 "<, > or <;> expected", &V);
+		if (*V.sym == comma)
+			 GetSym__(&V); 
+   } while (P_inset(*V.sym, mysys)); /* 11 */
+  
+
     test__(P_expset(SET, 1 << ((long)semicolon)), "<;> expected", &V);
+    };
     GetSym__(&V);
     globval.Elem_nFam++;
     if (globval.Elem_nFam <= Elem_nFamMax) {
       WITH = &ElemFam[globval.Elem_nFam-1];
+      //Aperture_Alloc(&WITH->ElemF);
       WITH1 = &WITH->ElemF;
+      Aperture_Alloc(WITH1);
       memcpy(WITH1->PName, ElementName, sizeof(partsName));
-      WITH1->PL = *V.rnum;
-      WITH1->Pkind = PartsKind(Aperture);
-      Aperture_Alloc(&WITH->ElemF);
-      AssignAPER(globval.Elem_nFam, &V);      
+      WITH1->Pkind = Aperture;
+      // WITH1->PL = QK; /* Length */ // need to track through a drift
+      AssignAPER(globval.Elem_nFam, &V);
+      CheckAperture(globval.Elem_nFam, &V); 
+
     } else {
       printf("Elem_nFamMax exceeded: %ld(%ld)\n",
 	     globval.Elem_nFam, (long)Elem_nFamMax);
       exit_(1);
     }
     break;
+    
 /**************************************************************************
       Marker
     ***************************************************************************
-- 
GitLab