/********************************************************/
/* ykj.c This file contains  functions for coordinate	*/
/* conversion YKJ coordinates.							*/
/* Alexander Kolesnikov.				*/
/********************************************************/


#include "ykj.h"


#include "geo_def.h"
#include "mercator.h"


//static const int nCorrection=1;


/*------------------------------------------------------*/


void YKJ_ForwardTransform(TCoordinate *Coordinate, TProjected *Projected,
						  int nCorrection)
{
  double lat, lon;
  double scale_fact;
  double northVal, eastVal;
  double lon_origin;
  double north_false, east_false;
  TDatum datum;
  TCoordinate coord_1;
  double north, east, DX, DY, a, b;
  TProjected proj_tmp;

  lat = Coordinate->Latitude;
  lon = Coordinate->Longitude;
  scale_fact = 1;

  lon_origin  = 27.0;
  north_false = 0;         /* False Northing */
  east_false  = 3500000;   /* False Easting */

  /*Make sure the longitude is between -180.00 .. 179.9*/
  lon = (lon+180) - (int)((lon+180)/360)*360-180;

  datum = International;
  if (nCorrection == 0) datum = International; /* KKJ(YKJ) uses Hayford 1910 ellipsoid */
  if (nCorrection == 1) datum = WGS84;

  coord_1.Latitude  = lat;
  coord_1.Longitude = lon;
  coord_1.Datum     = datum;

    ForwardConversion(&coord_1, lon_origin, &proj_tmp);
    eastVal  = east_false  + scale_fact * proj_tmp.Easting;
    northVal = north_false + scale_fact * proj_tmp.Northing;
   
  north=0; east=0;

 
  /* 0. No correction */
  if(nCorrection == 0)
  {
    north = northVal;
    east  = eastVal;
  }

  /* 4.1. 2D Rotation */
  if(nCorrection == 1)
  {
    DX = 132.858;      DY = 132.490;
    a  = 1.0000021251; b  = 0.0000044173;

    north = DX + a*northVal - b*eastVal;
    east  = DY + b*northVal + a*eastVal;
  }

  Projected->Northing   = north;
  Projected->Easting    = east;
  Projected->Datum      = International;
  Projected->ZoneNumber = 0;
}

/*------------------------------------------------------*/

void YKJ_InverseTransform(TProjected *Projected, TCoordinate *Coordinate,
						  int nCorrection)
{
  double north, east;
  double scale_fact;
  double lat_origin, lon_origin;
  double north_false, east_false;
  TDatum datum;
  double A,B,n1,e1,a1,b1,c,dN,dE;
  TProjected proj_tmp;

  /* Input parameters */
  north = Projected->Northing;   /* Northing */
  east  = Projected->Easting;    /* Easting */
  scale_fact = 1.0;             /* scale factor */

  lat_origin = 0.0;
  lon_origin = 27.0;
  east_false = 3500000.0;
  north_false= 0.0;

  /* 2D Rotation */
  datum = International;
  if(nCorrection == 0) datum = International;

  if(nCorrection == 1)
  {
    datum = WGS84;
    a1 = 1.0000021251; dN = 132.858;
    b1 = 0.0000044173; dE = 132.490;

    c  = a1*a1 + b1*b1;
    B  = b1/c;
    A  = a1/c;
    n1 =  A*(north - dN) + B*(east - dE);
    e1 = -B*(north - dN) + A*(east - dE);
    north = n1;
    east  = e1;
  }

  north = north - north_false;
  east  = east - east_false;

  proj_tmp.Northing   = north;
  proj_tmp.Easting    = east;
  proj_tmp.Datum      = datum;
  proj_tmp.ZoneNumber = 0;

  InverseConversion(&proj_tmp, lat_origin, lon_origin, scale_fact, Coordinate);
}

/*------------------------------------------------------*/

int YKJ_calcZoneNum(TCoordinate *Coordinate)
{
  int zone;

  zone = (int)((Coordinate->Longitude - 19.5)/3) + 1;
  if (zone < 1 || zone > 4) zone = -1;
  return zone;
}