/********************************************************/
/* utm.c This file contains  functions for coordinate	*/
/* conversion UTM coordinates.							*/
/* Alexander Kolesnikov.				*/
/********************************************************/



#include "utm.h"


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


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

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

static const double scale_fact=0.9996;

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

void UTM_ForwardTransform(TCoordinate *Coordinate, TProjected *Projected)
{
  double lat, lon;
  double northVal, eastVal;
  double lon_origin;
  double north_false, east_false;
  TCoordinate tmpCoord;
  TProjected proj_tmp;
  int intZoneNum; /* Zone: 1..60 */

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

  north_false = 0;      /* False Northing */
  east_false  = 500000; /* False Easting */

  if (lat < -80 || lat > 84)
  {
    printf("UTM: Imposible value of latitude\n");
    return;
  }

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

  tmpCoord.Latitude  = lat;
  tmpCoord.Longitude = lon;
  tmpCoord.Datum     = Coordinate->Datum;
  intZoneNum  = UTM_calcZoneNum(&tmpCoord);

  /* Set origin of the Zone */
  lon_origin  = intZoneNum*6 - 183;

  if(lat < 0) north_false = 10000000;   /*10000000 meter offset for southern hemisphere*/

    ForwardConversion(&tmpCoord, lon_origin, &proj_tmp);
    eastVal  = east_false  + scale_fact * proj_tmp.Easting;
    northVal = north_false + scale_fact * proj_tmp.Northing;

  Projected->Northing  = northVal;
  Projected->Easting    = eastVal;
  Projected->Datum      = proj_tmp.Datum;
  Projected->ZoneNumber = intZoneNum;
}

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

void UTM_InverseTransform(TProjected *Projected, TCoordinate *Coordinate)
{
  double lat_origin, lon_origin;
  double north_false, east_false;
  TProjected proj_tmp;
  int intZoneNum; /*Zone: 1..60 */
  double east, north;

  intZoneNum = Projected->ZoneNumber;
  east       = Projected->Easting;
  north      = Projected->Northing;

  lat_origin  = 0.0;
  east_false  = 500000.0;

  /* Set Zone parameters */

  int nz = (intZoneNum >= 0) ? intZoneNum : -intZoneNum;
  lon_origin  = ((6 * nz) - 183);


  if (intZoneNum < 0) north_false = 10000000.0;
  else north_false = 0.0;

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

  proj_tmp.Northing   = north;
  proj_tmp.Easting    = east;
  proj_tmp.Datum      = Projected->Datum;
  proj_tmp.ZoneNumber = intZoneNum;

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

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

int UTM_calcZoneNum(TCoordinate *Coordinate)
{
  int zone;
  double lat, lon;

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

  zone = (int)((lon + 180)/6) + 1;

  if (lat >= 56.0 && lat < 64.0 && lon >= 3.0 && lon < 12.0) zone = 32;

  /* Special zones for Svalbard */
  if (lat >= 72.0 && lat < 84.0)
  {
    if (lon >= 0.0 && lon <  9.0) zone = 31;
    else if (lon >= 9.0 && lon < 21.0) zone = 33;
         else if (lon >= 21.0 && lon < 33.0) zone = 35;
	      else if (lon >= 33.0 && lon < 42.0) zone = 37;
  }
  return zone;
}