AnalyticCoordinateTransformator.cc

Go to the documentation of this file.

#include <utl/AnalyticCoordinateTransformator.h>
#include <cmath>

using namespace std;
using namespace utl;


double
AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(const unsigned int i, const unsigned int j,
                                                                     const double r, const double theta, const double phi,
                                                                     const double e00, const double e11, const double e22,
                                                                     const double e01, const double e02, const double e12)
{
  const double t1 = cos(phi);
  const double t2 = t1 * t1;
  const double t3 = e00 * t2;
  const double t4 = cos(theta);
  const double t5 = t4 * t4;
  const double t6 = e00 * t5;
  const double t7 = t6 * t2;
  const double t8 = sin(theta);
  const double t9 = t8 * t2;
  const double t10 = r * t4;
  const double t13 = 0.2e1 * t9 * t10 * e01;
  const double t14 = t1 * r;
  const double t15 = sin(phi);
  const double t16 = t15 * e02;
  const double t18 = 0.2e1 * t14 * t16;
  const double t21 = 0.2e1 * t14 * t16 * t5;
  const double t22 = r * r;
  const double t23 = t22 * t5;
  const double t25 = t23 * t2 * e11;
  const double t26 = t22 * t4;
  const double t28 = t8 * t15;
  const double t31 = 0.2e1 * t26 * t1 * t28 * e12;
  const double t32 = t22 * e22;
  const double t33 = t22 * t2;
  const double t34 = t33 * e22;
  const double t37 = t33 * e22 * t5;
  const double t38 = t3 - t7 + t13 - t18 + t21 + t25 - t31 + t32 - t34 - t32 * t5 + t37;
  const double t39 = t1 * t15;
  const double t42 = t8 * t1;
  const double t44 = t4 * t15;
  const double t48 = t2 * r;
  const double t59 = r * e02;
  const double t63 = t22 * t15;
  const double t64 = t1 * e22;
  const double t68 = t39 * e00 - t39 * t6 + 0.2e1 * t42 * r * t44 * e01 + 0.2e1 * t48 * e02 - 0.2e1 * t48 * e02 * t5 + t23 * t39 * e11 + 0.2e1 * t26 * t9 * e12 - t59 + t59 * t5 - t26 * t8 * e12 - t63 * t64 + t63 * t64 * t5;
  const double t69 = t4 * e00;
  const double t72 = r * t5;
  const double t78 = r * t8;
  const double t82 = e12 * t5;
  const double t84 = t42 * t69 - t14 * e01 + 0.2e1 * t72 * t1 * e01 - t26 * t42 * e11 - t78 * t44 * e02 + t63 * e12 - t63 * t82;
  const double t87 = 0.2e1 * t10 * t8 * e01;
  const double t88 = t22 * e11;
  const double t89 = t88 * t5;
  const double t90 = e00 - t3 - t6 + t7 + t87 - t13 + t18 - t21 + t89 - t25 + t31 + t34 - t37;
  const double t102 = t22 * t1;
  const double t105 = t28 * t69 - r * t15 * e01 + 0.2e1 * t72 * t15 * e01 - t26 * t28 * e11 + t78 * t1 * t4 * e02 - t102 * e12 + t102 * t82;
  const double cg[][3] = {
    { t38, t68, t84 },
    { t68, t90, t105 },
    { t84, t105, t6 - t87 + t88 - t89 }
  };
  return cg[i][j];
}


double
AnalyticCoordinateTransformator::GetSphericalCovarianceFromCartesian(const unsigned int i, const unsigned int j,
                                                                     double x2, double y2, const double z2,
                                                                     const double e00, const double e11, const double e22,
                                                                     const double e01, const double e02, const double e12)
{
  if (fabs(x2) < 1e-10 && fabs(y2) < 1e-10) {
    x2 = 1e-5*z2;
    y2 = 1e-5*z2;
  }

  const double t1 = x2 *  x2;
  const double t2 = e00 *  t1;
  const double t5 = 2 * y2 * x2 * e01;
  const double t6 = x2 * z2;
  const double t9 = y2 * y2;
  const double t10 = e11 * t9;
  const double t11 = y2 * z2;
  const double t14 = z2 * z2;
  const double t17 = t1 + t9 + t14;
  const double t18 = 1 / t17;
  const double t20 = t1 * x2;
  const double t22 = z2 * t1;
  const double t29 = y2 * e01;
  const double t32 = x2 * t14;
  const double t34 = z2 * t9;
  const double t36 = t9 * y2;
  const double t42 = t20 * e02 - e00 * t22 + e22 * t22 + e12 * y2 * t1 + t9 * e02 * x2 - 2 * t29 * t6 - e02 * t32 - e11 * t34 + e12 * t36 + t9 * e22 * z2 - e12 * y2 * t14;
  const double t43 = t1 + t9;
  const double t45 = sqrt( (t18 * t43));
  const double t46 = 0.1e1 / t45;
  const double t48 = t17 * t17;
  const double t49 = 1 / t48;
  const double t50 =  t49 * t46 *  t42;
  const double t61 = sqrt( t17);
  const double t64 = 1 / t43;
  const double t65 =  t64 / t61 *  (-y2 * e00 * x2 + t1 * e01 - e01 * t9 + y2 * e11 * x2 - e02 * t11 + z2 * e12 * x2);
  const double t66 = t1 * t1;
  const double t84 = t9 * t9;
  const double t110 =  t64 * t46 / t61 /  t17 *  (e00 * x2 * t11 - t1 * e01 * z2 + e01 * t34 - e11 * x2 * t11 - t1 * e02 * y2 - e02 * t36 + t20 * e12 + t9 * x2 * e12);
  const double t114 = t43 * t43;
  const double cg[][3] = {
    { t18 * (t2 + t5 + 2 * e02 * t6 + t10 + 2 * e12 * t11 + t14 * e22), -t50, t65 },
    { -t50, t49 * t64 * (t66 * e22 - 2 * t20 * e02 * z2 - 2 * y2 * e12 * t22 + 2 * t9 * e22 * t1 + t14 * t2 - 2 * t9 * e02 * t6 + 2 * t32 * t29 + t14 * t10 + t84 * e22 - 2 * e12 * t36 * z2), -t110 },
    { t65, -t110, 1 / t114 * (e00 * t9 - t5 + e11 * t1) }
  };
  return cg[i][j];
}


double
AnalyticCoordinateTransformator::GetRadiusFromCartesianCoordinates(const double x, const double y, const double z)
{
  return sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2));
}


double
AnalyticCoordinateTransformator::GetAzimuthFromCartesianCoordinates(const double x, const double y)
{
  double azimuth = atan2(y, x);
  if (azimuth < 0) {
    azimuth += 2*M_PI;
  }
  return azimuth;
}


double
AnalyticCoordinateTransformator::GetZenithFromCartesianCoordinates(const double x, const double y, const double z)
{
  return acos(z / GetRadiusFromCartesianCoordinates(x,y,z));
}


double
AnalyticCoordinateTransformator::GetXFromSphericalCoordinates(const double r, const double zenith, const double azimuth)
{
  return r*cos(azimuth)*sin(zenith);
}


double
AnalyticCoordinateTransformator::GetYFromSphericalCoordinates(const double r, const double zenith, const double azimuth)
{
  return r*sin(azimuth)*sin(zenith);
}


double
AnalyticCoordinateTransformator::GetZFromSphericalCoordinates(const double r, const double zenith)
{
  return r*cos(zenith);
}