Changeset 3430 in Sophya

Timestamp:

Dec 11, 2007, 10:32:45 AM (18 years ago)

Author:

ansari

Message:

Ajout commande errorellipse (fait pour these de S. Bargot) a sopiamodule.cc - Reza 11/12/2007

File:

• trunk/SophyaPI/ProgPI/sopiamodule.cc (modified) (6 diffs)

trunk/SophyaPI/ProgPI/sopiamodule.cc

@@ -32 +32 @@
 }
 
+// --- fonctions pour faire FFT
 void SophyaFFT_forw(string& nom, string& nomout, string dopt);
 void SophyaFFT_back(string& nom, string& nomout, string dopt, bool forcez8=false);
@@ -37 +38 @@
 void SophyaFFT_filter(string& nom, string& filt, string& outvec, string dopt, bool fgvec=false);
 
+//--- fonctions pour trace d'ellipse d'erreur
+int  diag_ellipse(double s1, double s2, double c, 
+                  double &a, double &b, double& tet, double& phi, int prt=0);
+
+void ellipse2vec(double a, double b, double tet, Vector& evx, Vector& evy, 
+                 double scale=1., int nbpoints=101);
+
+//--- Declaration du CmdExecutor du module
 class sopiamoduleExecutor : public CmdExecutor {
 public:
@@ -44 +53 @@
 };
 
-//  Classe pour trace de grille en projection spherique Mollweide 
+//---  Classe pour trace de grille en projection spherique Mollweide 
 class MollweideGridDrawer : public PIDrawer {
 public:
@@ -188 +197 @@
   SetLimits(x0-xlarg/2., x0+xlarg/2., y0-ylarg/2., y0+ylarg/2.);
 }
-
+// Cette fonction calcule les parametres de l'ellipse des erreurs
+// a,b 1/2 grand, petit axe ; tet : Angle du grand axe / Ox 
+// A partir de s1 = Sigma_1^2 , s2 = Sigma_2^2 , c = Covar_1_2 = Sigma_1 x Sigma_2 x Rho 
+
+//---------------------------------------------------------------------------
+//--- Fonctions de trace d'ellipses d'erreur 
+//---------------------------------------------------------------------------
+
+/* Nouvelle-Fonction */
+int  diag_ellipse(double s1, double s2, double c, 
+                  double &a, double &b, double& tet, double& phi, int prt)
+{
+  int rc = 0;
+  // On calcule les valeurs propres de la matrice 
+  // | s1  c  |
+  // | c   s2 |
+  a = b = 0.;  // Demi axes de l'ellipse
+  tet = 0.;     // Cos theta de l'angle de l'axe de l'ellipse
+  phi = 0.;
+  double det = (s1-s2)*(s1-s2) + 4*c*c;
+  if (det < 0.) {
+    if (prt > 0) 
+      cout << "ERREUR : diag_ellipse, ERREUR det = " << det << " <0. s1=" 
+           << s1 << " s2=" << s2 << " c=" << c << endl;
+    return 99;
+  }
+  // les deux valeurs propres l1, l2
+  double l1 = 0.5*((s1+s2)+sqrt(det));
+  double l2 = 0.5*((s1+s2)-sqrt(det));
+  if (l1 < 0.) {
+    if (prt > 1) 
+      cout << "diag_ellipse/Warning l1 < 0. -> a = -sqrt(l1)" << endl; 
+    a = -sqrt(-l1);
+    rc += 2;
+  }
+  else   a = sqrt(l1);
+  if (l2 < 0.) {
+    if (prt > 1) 
+      cout << "diag_ellipse/Warning l2 < 0. -> b = -sqrt(l2)" << endl; 
+    b = -sqrt(-l2);
+    rc += 4;
+  }
+  else b = sqrt(l2);
+  
+  if (fabs(c) > 1.e-9)  tet = atan2((l1-s1), c);
+  else {
+    if (prt > 2) 
+      cout << "diag_ellipse/Warning c (Correlation) = 0. -> Tet=0 " << endl; 
+  }
+  // Calcul de tan(2 phi)
+  double t2phi = 0.;
+  if (fabs(s1-s2) > 1.e-9) {
+    t2phi = 2.*c/(s1-s2);
+    phi = atan(t2phi)/2.;
+  }
+  else {
+    if (prt > 3) 
+      cout << "diag_ellipse/Warning fabs(s1-s2)=" << fabs(s1-s2)
+           << " -> Phi=Pi/4 ... " << endl; 
+    phi = 0.25*M_PI;
+  }
+  if (prt > 0) 
+    cout << "diag_ellipse/Info: s1=" << s1 << " s2=" << s2 
+         << " c=" << c << " ---> l1=" << l1 << " l2=" << l2 
+         << "\n a= " << a << " b=" << b << " Theta= " << tet 
+         << "( " << tet*2/M_PI << " xPi/2)  Phi= " << phi 
+         << "( " << phi*2/M_PI << " x Pi/2) " << endl;
+  return rc;
+}
+
+/* Nouvelle-Fonction */
+void ellipse2vec(double a, double b, double tet, Vector& evx, Vector& evy, 
+                 double scale, int nbpoints)
+{
+  a *= scale;
+  b *= scale;
+  
+  int npt = nbpoints;
+  double pas = M_PI*2./(npt-1);
+  evx.SetSize(npt+1);
+  evy.SetSize(npt+1);
+  for(int i=0; i<=npt; i++) {
+    double t = i*pas;
+    double xr = a*cos(t);
+    double yr = b*sin(t);
+    double ar = atan2(yr, xr);
+    double r = sqrt(xr*xr+yr*yr);
+    evx(i) = r*cos(ar+tet);
+    evy(i) = r*sin(ar+tet);
+  }
+  return;
+}
+
+//---------------------------------------------------------------------------
+//--- Classe sopiamoduleExecutor
+//---------------------------------------------------------------------------
 /* --Methode-- */
 sopiamoduleExecutor::sopiamoduleExecutor()
@@ -251 +355 @@
 usage = "Set default value for Molleweide projection of spherical maps (outside maps)";
 usage += "\n Usage: setprjmoldefval OutOfMapValue ";
+mpiac->RegisterCommand(kw, usage, this, hgrp);
+
+//--- Trace d'ellipse d'erreur
+hgrp = "SophyaCmd";
+kw = "errorellipse";
+usage = " Error ellipse plotter \n";
+usage += "Usage: errellipse [-fill] xc yc sx2 sy2 cov [scale=1] [gr_opt] [npt=128] \n"; 
+usage = "   Adds an error ellipse, specified by its center (xc,yc) \n";
+usage += "  and error matrix (covariance) parameters  SigX^2,SigY^2,Covar \n" ;
+usage += "  A global scaling parameters scale , graphic attributes and \n";
+usage += "   nomber of points can optionnaly be specified\n";
+usage += "  if the -fill flag specified, plots a filled ellipse \n" ;
 mpiac->RegisterCommand(kw, usage, this, hgrp);
 
@@ -358 +474 @@
    wrsid = omg.GetImgApp()->DispScDrawer(mollgrid, name, dopt); 
  }
+else if (kw == "errorellipse") {
+  if ((tokens.size()<5) || ((tokens.size()>0)&&(tokens[0]=="-fill")&&(tokens.size()<6))) {
+    cout << "Usage: errorellipse [-fill] xc yc sx2 sy2 cov [scale=1 gropt npt=128]" << endl;
+    return(1);
+  }
+  bool fgfill = false;
+  unsigned int offtok = 0;
+  if (tokens[0] == "-fill") {
+    offtok = 1;  fgfill = true;
+  }
+  double xc = atof(tokens[0+offtok].c_str());
+  double yc = atof(tokens[1+offtok].c_str());
+  double sx = atof(tokens[2+offtok].c_str());
+  double sy = atof(tokens[3+offtok].c_str());
+  double cov = atof(tokens[4+offtok].c_str());
+  double scale = 1.;
+  if (tokens.size() > (5+offtok)) scale = atof(tokens[5+offtok].c_str());
+  string gropt;
+  if (tokens.size() > (6+offtok)) gropt = tokens[6+offtok];
+  int npt = 128;
+  if (tokens.size() > (7+offtok)) npt = atoi(tokens[7+offtok].c_str());
+
+  double a,b,tet,phi;
+  diag_ellipse(sx, sy, cov, a, b, tet, phi, 1);
+  Vector evx, evy;
+  ellipse2vec(a,b,tet,evx, evy, scale, npt+1);
+  vector<double> vx, vy;
+  for(int k=0; k<evx.Size(); k++) {
+    vx.push_back(evx(k)+xc);
+    vy.push_back(evy(k)+yc);
+  }
+  NamedObjMgr omg;
+  omg.GetImgApp()->AddPoly(vx, vy, gropt, fgfill, false);
+}
 
 return(0);