Коническая равноугольная

//---------------------------------------------------------------------
// Инициализировать константы для пересчета координат
//---------------------------------------------------------------------
void TTranslate::InitConicalOrthomorphic()
{
  FalseEasting = 8000000;

  Sign = 1;

  if ((FirstMainParallel < 0) ||
      ((FirstMainParallel == 0) && (SecondMainParallel < 0)))
    {
      Sign = -1;
      FirstMainParallel  = -FirstMainParallel;
      SecondMainParallel = -SecondMainParallel;
    }

  double   a1 = M_PI;
  double kci1 = asin(ExscMerid * sin(FirstMainParallel));
  double kci2 = asin(ExscMerid * sin(SecondMainParallel));
  double   u1 = tan(a1/4 + FirstMainParallel/2)/pow(tan(a1/4+kci1/2),ExscMerid);
  double   u2 = tan(a1/4 + SecondMainParallel/2)/pow(tan(a1/4+kci2/2),ExscMerid);
  double   r1 = BigAxis * cos(FirstMainParallel)/sqrt(1-pow(ExscMerid*sin(FirstMainParallel),2));
  double   r2 = BigAxis * cos(SecondMainParallel)/sqrt(1-pow(ExscMerid*sin(SecondMainParallel),2));
           va = (log10(r1)-log10(r2))/(log10(u2)-log10(u1));

  if (va == 0)
    va = 1;

  vd = r1 * pow(u1, va)/va;

  // Перевод геодезических координат левого нижнего угла в
  // прямоугольные координаты

  kci1  = asin(ExscMerid * sin(MainPointParallel));
  u1    = tan(a1/4 + MainPointParallel/2)/pow(tan(a1/4+kci1/2),ExscMerid);

  vb = vd/pow(u1, va);

  vc = M_PI / 2.0 / va;
}

//---------------------------------------------------------------------
// Вычисление геодезических координат точки B,L(равноугольной конической
// проекции) по координатам точки X,Y
// Входные данные:
// x,y - координаты точки в метрах
// Выходные данные:
// b,l - широта и долгота точки в радианах
//---------------------------------------------------------------------
void _fastcall TTranslate::XY2BL_ConicalOrthomorphic(double x, double y, double& b, double& l)
{
 // Перевод координат точки из равноугольной проекции в
 // геодезические координаты
 // Ищем долготу

 double delta = 0.000000001;      // 1E-9

 y =  y - FalseEasting;
 x = (x - FalseNorthing) * Sign;

 if (fabs(y) > delta)
   {
     double dqx = vb - x;
     if (fabs(dqx) < delta)
       dqx = delta;

     if (dqx < 0)                  // Переход через Полюс
       {
         double zero = vc;     // atan(y/delta)/Dlog;
         if (y < 0)
           zero = -zero;

         l = zero + zero - atan(y/(-dqx))/va;
       }
     else
       {
         l = atan(y/dqx)/va;
       }
   }
 else
   {
     l=0.0;
     if (y > 0)
       y = delta;
     else
       y = -delta;
   }

 double sigma = va * l;

 l += AxisMeridian;

 // Ищем широту
 double d = -M_PI/4;
 double h =  M_PI/2-0.0001;

// h=0; d=-M_PI/2+0.0001;
 double rab1 = sin(sigma) * vd / y;
 if (fabs(rab1) <= delta)
   {
     b = M_PI/2.;
     return;
   }

 double fd, fh, c, fc;

 // Заглушка для левых значений
 if (rab1 < 0)
   {
     rab1 = 1;
     fd   = 0;
     fh   = 0;
   }
 else
   {
     rab1 = pow(rab1,1./ va);

     fd   = sqrt((1+sin(d))*pow(1-ExscMerid*sin(d),ExscMerid)/((1-sin(d))*pow(1+ExscMerid*sin(d),ExscMerid))) -
            rab1;
     fh   = sqrt((1+sin(h))*pow(1-ExscMerid*sin(h),ExscMerid)/((1-sin(h))*pow(1+ExscMerid*sin(h),ExscMerid))) -
            rab1;
   }

 if (fd == 0)
   {
     b = d;
     b = b * Sign;
     return;
   }

 if (fh == 0)
   {
     b = h;
     b = b * Sign;
     return;
   }

 while (h-d > delta)
 {
   if (fd*fh < 0)
     {
       c  = (h-d)/2;
       c +=d;

       fc = sqrt((1+sin( c ))*pow(1-ExscMerid*sin( c ),ExscMerid)/((1-sin( c ))*pow(1+ExscMerid*sin( c ),ExscMerid))) - rab1;

       if (fabs(fc) < delta)
         {
           b = c;
           b = b * Sign; //24/02/04
           return;
         }

       if (fd*fc < 0)
         h = c;
       else
         {
           d  = c;
           fd = fc;
         }
     }
   else
     {
       b = h;
       b = b * Sign;
       return;
     }

   b = c;
 }

 b = b * Sign;
}

//---------------------------------------------------------------------
// Вычисление координат X,Y точки (равноугольной
// конической проекции) по геодезическим координатам точки B,L
// Входные данные:
// b,l - геодезические координаты точки в радианах
// Выходные данные:
// x,y - координаты точки в метрах
//---------------------------------------------------------------------
void _fastcall TTranslate::BL2XY_ConicalOrthomorphic(double b, double l, double& x, double& y)
{
 double a1 = M_PI;

 // Перевод геодезических координат точки в прямоугольные координаты
 b = b * Sign;

 double kci1 = asin(ExscMerid * sin(b));
 double u1 = tan(a1/4 + b/2)/pow(tan(a1/4+kci1/2),ExscMerid);

 double p = pow(u1, va);                     // va -> Dlog

 if (p < DOUBLENULL)                         // 10**-6
   p = 1;

 p = vd/p;vd -> K

 double sigma = va * (l - AxisMeridian);

 y = p * sin(sigma);
 x = vb - p * cos(sigma);

 // Перевод координат точки в систему листа
 y = y + FalseEasting;
 x = x * Sign + FalseNorthing;
}

//--------------------------------------------------------------------
// Инициализация параметров эллипсоида
//--------------------------------------------------------------------
void TTranslate::InitEllipsoid(int ellipsoid, ELLIPSOIDPARAM * ellparam)
{
  if ((ellipsoid == USERELLIPSOID) && (ellparam != 0) &&
      (ellparam->SemiMajorAxis > 0) && (ellparam->InverseFlattening >= 0))
    {
       BigAxis   = ellparam->SemiMajorAxis;       
       Alfa      = ellparam->InverseFlattening;
    }
  else
    {
       if ((ellipsoid <= 0) || (ellipsoid > ELLIPSOIDCOUNT))
         ellipsoid = WGS_84;

       BigAxis   = Spheroids[ellipsoid].SemiMajorAxis;
       Alfa      = Spheroids[ellipsoid].InverseFlattening;
    }

  Ellipsoid = ellipsoid;

  AlfaTo1   = 1.0 - Alfa;

  E2        = 2.*Alfa - Alfa * Alfa;
  if (E2 > 0)
    ExscMerid = sqrt(E2);
  else
    ExscMerid = 0;
}