Changeset 341 in PSPA for Interface_Web/trunk/pspaWT/sources/controler/src/particleBeam.cc

Timestamp:

Feb 22, 2013, 10:07:46 AM (11 years ago)

Author:

lemeur

Message:

fin dessins esp. phase + noms de fichiers

File:

• Interface_Web/trunk/pspaWT/sources/controler/src/particleBeam.cc (modified) (1 diff)

Interface_Web/trunk/pspaWT/sources/controler/src/particleBeam.cc

@@ -426 +426 @@
 }
 
+
+
+void particleBeam::donneesDessinEllipse(vector<double>& xcor, vector<double>& ycor, unsigned indexAbs, unsigned indexOrd) {
+  int k;
+  double x,y;
+
+  if ( !momentRepresentationOk_ ) return;
+
+  if ( indexAbs > 5 || indexOrd > 5 ) return;
+
+  xcor.clear();
+  ycor.clear();
+  // les index sont dans l'ordre x,y,z,xp,yp, de/E
+  // on traduit en TRANSPORT
+  if ( indexAbs == 1 ) indexAbs = 2;  // y
+  if ( indexAbs == 2 ) indexAbs = 4;  // z -> l
+  if ( indexAbs == 3 ) indexAbs = 1;  // xp
+  if ( indexAbs == 4 ) indexAbs = 3;  // yp
+
+  if ( indexOrd == 1 ) indexOrd = 2;  // y
+  if ( indexOrd == 2 ) indexOrd = 4;  // z -> l
+  if ( indexOrd == 3 ) indexOrd = 1;  // xp
+  if ( indexOrd == 4 ) indexOrd = 3;  // yp
+
+  cout << " index x" << indexAbs << " index y " << indexOrd << endl;
+
+  double xm = ( rij_.getMatrix().at(indexAbs) ).at(indexAbs);
+  double ym = ( rij_.getMatrix().at(indexOrd) ).at(indexOrd);
+  double r;
+  if ( indexOrd > indexAbs ) {
+    r  = ( rij_.getMatrix().at(indexOrd) ).at(indexAbs);
+  } else {
+    r  = ( rij_.getMatrix().at(indexAbs) ).at(indexOrd);
+  }
+
+  cout <<  " racs11= " << xm << " racs22= " << ym << " r12= " << r << endl;
+
+
+  int nbintv = 50;
+  if ( xm == 0.0 ) return;
+  double pas = 2.0 * xm / nbintv;
+
+  //  cout << " r= " << r << endl;
+  double rac = (1 - r*r);
+  if ( rac > 0.0 ) 
+    {
+      cout << " cas rac > " << endl;
+      rac = sqrt(1 - r*r);
+      double alpha = -r / rac;
+      double beta = xm / ( ym * rac);
+      //  double gamma = ym / ( xm * rac );
+      double epsil = xm * ym * rac;
+      double fac1 = -1.0 / ( beta * beta);
+      double fac2 = epsil/beta;
+      double fac3 = -alpha/beta;
+      double aux;
+      for ( k=0; k < nbintv; k++)
+        {
+          x = -xm + k*pas;
+          aux = fac1 * x * x + fac2;
+          //     cout << " aux2= " << aux << endl;
+          if ( aux <= 0.0 )
+            {
+              aux = 0.0;
+            }
+          else aux = sqrt(aux);
+      
+          //        y = fac3*x;
+          y = fac3*x + aux;
+          xcor.push_back(x);
+          ycor.push_back(y);
+        }
+
+      for ( k=0; k <= nbintv; k++)
+        {
+          x = xm - k*pas;
+          aux =  fac1 * x * x + fac2;
+          if ( aux <= 0.0 ) 
+            {
+              aux = 0.0;
+            }
+          else aux = sqrt(aux);
+          //   y = fac3*x;
+          y = fac3*x - aux;
+          xcor.push_back(x);
+          ycor.push_back(y);
+        }
+    }
+  else
+    // cas degenere
+    {
+      cout << " cas degenere " << endl;
+      double fac = ym/xm;
+      for ( k=0; k < nbintv; k++)
+        {
+          x = -xm + k*pas;
+          y = fac*x;
+          xcor.push_back(x);
+          ycor.push_back(y);
+        }
+        
+    }
+}
+
 void particleBeam::histogramme(vector<double>&xcor,vector<int>& hist,int& cnts,double out[3])
 {