/********************************************************/
/* kkj.c This file contains  functions for coordinate	*/
/* conversion KKJ coordinates.							*/
/* Alexander Kolesnikov.				*/
/********************************************************/


#include "kkj.h"


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

#include <stdio.h>
#include <math.h>


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

//static const int nCorrection = 1; 

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

void KKJ_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;
  int         intZoneNum, nZ;

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

  north_false = 0;     /* False Northing */

  /*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;
  if (nCorrection == 1) datum = WGS84;

  if (nCorrection == 2) 
  {
	  datum = International;       // datum_ED50;
      // !!!! coord_1 = Coordinate.WGS84_to_ED50(coord_0);
  }
 
  
  intZoneNum = KKJ_calcZoneNum(Coordinate);

  /* Calc false easting */
  /* 1 500 000; 2 500 000; 3 500 000; 4 500 000 */
  east_false = 1500000 + (intZoneNum-1)*1000000;

  /* Set origin of the Zone */
  lon_origin = 21 + (intZoneNum-1)*3;

  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.  Correction: 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;
  }
 
  // 4.2. 2D+3D  Rotations;
  if(nCorrection == 2)
  {
     DX = -61.5805;      DY = 95.6691;
     a = 1.0000007515;   b = 0.0000043933;
     
	 nZ = KKJ_calcZoneNum(Coordinate);

     north = DX + a*northVal - b*(eastVal-nZ*1000000.0);
     east  = DY + b*northVal + a*(eastVal-nZ*1000000.0) + nZ*1000000.0;
  }

 /*******************************************/

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

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

void KKJ_InverseTransform(TProjected *Projected, TCoordinate *Coordinate, int nCorrection)
{
  int        intZoneNum, nZ; /*Zone: 1..7 */
  double     north, east;
  double     scale_fact;
  double     lat_origin, lon_origin;
  double     north_false, east_false;
  double     a, b, A,B,n1,e1,a1,b1,c,dN,dE;
  TDatum     datum;
  TProjected proj_tmp;
  
  north      = Projected->Northing;   /* Northing */
  east       = Projected->Easting;    /* Easting */
  intZoneNum = Projected->ZoneNumber;
  scale_fact = 1.0;                  /* scale factor */

  lat_origin = 0.0;
  north_false= 0.0;

  /* Set Zone parameters */

  int nz =   (intZoneNum >= 0) ? intZoneNum  : -intZoneNum ;
//  lon_origin = 21 + 3*(abs(intZoneNum)-1);
  lon_origin = 21 + 3*(nz-1);

  east_false = 1500000 + (intZoneNum-1)*1000000;

  /* 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;
  }

  if(nCorrection == 2)
  {
      a = 1.0000007515;   dN = -61.5805;
      b = 0.0000043933;   dE = 95.6691;
      datum = International;
      nZ    = (int)(east/1000000.0);

	  c  = a*a + b*b;
      B  = b/c;
      A  = a/c;
      n1 =  A*(north - dN) +  B*(east - dE - nZ*1000000);
      e1 = -B*(north - dN) +  A*(east - dE - nZ*1000000) + nZ*1000000;

      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);


  /* !!!!!!!!!
  if(nCorrection == 2) 
  {
      coord_dst  = Coordinate.ED50_to_WGS84(coord_tmp);
  }
  else
  {
      coord_dst = coord_tmp;
  }

*/



}

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

int KKJ_calcZoneNum(TCoordinate *Coordinate)
{
  int zone;

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