/*******************************************************************************
/ Código escrito para generar un binario que obtenga las longitudes y las
/ latitudes de las esquinas de las capas de los mapas para hacer unas pruebas
/ en Android.
/
/ Marcos Molina Cano - 10/02/2016
*******************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

int orthographic_inv ( double *lambda, double *phi, double x, double y, double lambda0, double phi0) {//, double *vaux ) {
  double r, c, xx, yy;
  r = sqrt ( x * x + y * y );
  if ( r >= 1.0 )
    return 0;
  c = asin ( r );
  if ( r > 0.0 ) {
      *phi = asin ( cos ( c ) * sin ( phi0 ) + y * sin ( c ) * cos ( phi0 ) / r );
      *lambda = lambda0 +  atan ( x * sin ( c ) / ( r * cos ( phi0 ) * cos ( c ) - y * sin ( phi0 ) * sin ( c ) ) );
    }
  else {
    *phi = phi0;
    *lambda = lambda0;
    return 1;
  }
}

int mercator_inv ( double *lambda, double *phi, double x, double y, double lambda0, double phi0) { // double *vaux )
  *lambda = x  + lambda0;
  *phi = atan(sinh( y )); // + phi0;

  if ( *lambda >= M_PI ) *lambda -= M_PI * 2;
  if ( *lambda < ( -M_PI ) ) *lambda += M_PI * 2;
//  *vaux = phi0; // para evitar warnings
  return 1;
}

int main(int argc, char const *argv[]) {
  struct map {
    char CPROY; /*!< caracter que codifica la proyeccion. O=ortográfica, G=geostacionaria, L=lambert */
    double XLAT0; /*!< latitud central de la proyeccion, radianes */
    double XLON0; /*!< longitud central de la proyeccion, radianes */
    double XLAT; /*!< latitud central del mapa, radianes */
    double XLON; /*!< longitud central del mapa, radianes */
    double ESCALA; /*!< escala del mapa */
    double FACTORX; /*!< factor de escala X */
    double FACTORY; /*!< factor de escala Y */
    double Dx; /*!< desplazamiento centro grafico coor X respecto centro proy */
    double Dy; /*!< desplazamiento centro grafico coor Y respecto centro proy */
    int SIZEX; /*!< anchura en pixels del mapa */
    int SIZEY; /*!< altura en pixels del mapa */
    int X0; /*!< coordenada X del pixel central*/
    int Y0; /*!< coordenada Y del pixel central */
  } MAPA;

  // mapa_espana -------------> "penin"
  // mapa_europa -------------> "europ"
  // mapa_canarias -----------> "canar"
  // mapa_francia ------------> "france"
  // mapa_alemania -----------> "aleman"
  // mapa_inglaterra ---------> "ukn"
  // mapa_italia -------------> "italia"
  // mapa_portugal -----------> "portugal"
  // mapa_madeira ------------> "madeira"
  // mapa_chile_norte --------> "acn"
  // mapa_chile_sur ----------> "acs"
  // mapa_argentina_norte ----> "acn"
  // mapa_argentina_sur ------> "acs"
  // mapa_brasil_amazonia ----> "bra"
  // mapa_brasil_este --------> "bre"
  // mapa_brasil_sur ---------> "brs"
  // mapa_mexico_oeste -------> "mexw"
  // mapa_mexico_este --------> "mexe"
  // mapa_usa_oeste ----------> "usaw"
  // mapa_usa_este -----------> "usae"
  // mapa_honduras -----------> "amercen"
  // mapa_colombia -----------> "vencol"
  // mapa_canada_este --------> "cnde"
  // mapa_canada_oeste -------> "cndw"
  // mapa_rusia --------------> "rusia"
  // mapa_siberia_occidental -> "sibocc"
  // mapa_siberia_oriental ---> "sibori"
  // mapa_pb (paises bajos) --> "holand"
  // mapa_austria ------------> "austria"
  // mapa_oceania ------------> "oce"
  // mapa_africa -------------> "afr"
  // mapa_america_central ----> "amc"
  // mapa_america_sur --------> "ams"
  // mapa_asia ---------------> "asi"

  char *nombresZonas[34] = {
    "espana", "europa", "canarias", "francia", "alemania", "inglaterra",
    "italia", "portugal", "madeira", "chile_norte", "chile_sur", "argentina_norte",
    "argentina_sur", "brasil_amazonia", "brasil_este", "brasil_sur", "mexico_oeste", "mexico_este",
    "usa_oeste", "usa_este", "honduras", "colombia", "canada_este", "canada_oeste",
    "rusia", "siberia_occidental", "siberia_oriental", "pb", "austria", "oceania",
    "africa", "america_central", "america_sur", "asia"
  };

  char *zonas[34] = { "penin", "europ", "canar", "france", "aleman", "ukn",
                      "italia", "portugal", "madeira", "acn", "acs", "acn",
                      "acs", "bra", "bre", "brs", "mexw", "mexe",
                      "usaw", "usae", "amercen", "vencol", "cnde", "cndw",
                      "rusia", "sibocc", "sibori", "holand", "austria", "oce",
                      "afr", "amc", "ams", "asi" };

//  char *zonasMercator[31] = { "africa", "alaska", "alemania", "america_central", "america_sur", "argentina",
//                            "asia", "austria", "azores", "brasil", "canada", "canarias", "chile", "colombia",
//                            "eeuu", "espana", "europa", "francia", "hawaii", "honduras", "india", "inglaterra", "italia",
//                            "madeira", "mexico", "oceania", "paises_bajos", "portugal", "rusia", "siberia_occidental",
//                            "siberia_oriental" };

//  char *zonasMercator[4] = {"kazajistan", "nueva_zelanda", "rusia2", "marruecos"};
//  char *zonasMercator[2] = {"belgica", "dinamarca"};
  char *zonasMercator[1] = {"polonia"};

  double pi = 4.0 * atan( 1.0 );
  double x, y, lat, lon;
  int ix, iy, k;

  FILE *mtfFile, *destinyFile;
  char nombre[255];

  destinyFile = fopen( "/home/marcos/wrf_model/varios/lat-lon_mercator4.txt", "w" );
  if ( destinyFile == NULL) {
    printf( "*** Error *** --> No se ha podido crear destinyFile\n" );
  }

  for ( k = 0; k < 30; k++ ) {
    /* Nombre del archivo *.mtf con su ruta. */
//    strcpy( nombre, "/home/marcos/meteored/share/" );
//    strcat( nombre, zonas[k] );
//    strcat( nombre, "_lyr.mtf");
      strcpy( nombre, "/home/marcos/polonia/" );
      strcat( nombre, zonasMercator[k]);
      strcat( nombre, "_M.mtf");


    mtfFile = fopen( nombre, "r" );
    if ( mtfFile == NULL) {
      printf( "*** Error *** --> No se ha podido abrir %s. Saliendo...\n", nombre );
      return -1;
    }
    fscanf( mtfFile, "%c", &MAPA.CPROY );
    fscanf( mtfFile, "%lf %lf", &MAPA.XLAT0, &MAPA.XLON0 );
    fscanf( mtfFile, "%lf %lf", &MAPA.XLAT, &MAPA.XLON );
    fscanf( mtfFile, "%lf", &MAPA.ESCALA );
    fscanf( mtfFile, "%lf %lf", &MAPA.FACTORX, &MAPA.FACTORY );
    fscanf( mtfFile, "%lf %lf", &MAPA.Dx, &MAPA.Dy );
    fscanf( mtfFile, "%i %i", &MAPA.SIZEX, &MAPA.SIZEY );
    fscanf( mtfFile, "%i %i", &MAPA.X0, &MAPA.Y0 );

    // printf("%c\n", MAPA.CPROY );
    // printf("%lf %lf\n", MAPA.XLAT0, MAPA.XLON0 );
    // printf("%lf %lf\n", MAPA.XLAT, MAPA.XLON );
    // printf("%lf\n", MAPA.ESCALA );
    // printf("%lf %lf\n", MAPA.FACTORX, MAPA.FACTORY );
    // printf("%lf %lf\n", MAPA.Dx, MAPA.Dy );
    // printf("%i %i\n", MAPA.SIZEX, MAPA.SIZEY );
    // printf("%i %i\n", MAPA.X0, MAPA.Y0 );

    fclose( mtfFile );

    fprintf( destinyFile, "%s%s\n", "mapa_", zonasMercator[k] );

    /* Vértice inferior izquierdo. */
    ix = 0;
    iy = MAPA.SIZEY - 1;
    x = ( double ) ( ix - MAPA.X0 ) / MAPA.FACTORX + MAPA.Dx;
    y = ( double ) ( MAPA.Y0 - iy ) / MAPA.FACTORY + MAPA.Dy;
    // printf("x=%lf y=%lf\n", x, y);
    if ( MAPA.CPROY == 'O' )
      orthographic_inv( &lon, &lat, x, y, MAPA.XLON0, MAPA.XLAT0);
    else if ( MAPA.CPROY == 'M')
      mercator_inv( &lon, &lat, x, y, 0.0, 0.0 );
    else {
      printf("*** ERROR *** - No se reconoce la proyección elegida.");
      return -1;
    }
    if ( lon > 180.0 ) {
      lon -= 2.0 * pi;
    }
    printf("Vértice inferior izquierdo:\n" );
    fprintf( destinyFile, "%s", "Vértice inferior izquierdo:\n" );
    printf("  lon=%lf lat=%lf\n", lon*180.0/pi, lat*180/pi );
    fprintf( destinyFile, "  %s%lf %s%lf\n", "lon=", lon*180.0/pi, "lat=", lat*180/pi);

    /* Vértice inferior derecho. */
    ix = MAPA.SIZEX - 1;
    iy = MAPA.SIZEY - 1;
    x = ( double ) ( ix - MAPA.X0 ) / MAPA.FACTORX + MAPA.Dx;
    y = ( double ) ( MAPA.Y0 - iy ) / MAPA.FACTORY + MAPA.Dy;
    // printf("x=%lf y=%lf\n", x, y);
    if ( MAPA.CPROY == 'O' )
      orthographic_inv( &lon, &lat, x, y, MAPA.XLON0, MAPA.XLAT0);
    else if ( MAPA.CPROY == 'M')
      mercator_inv( &lon, &lat, x, y, 0.0, 0.0 );
    else {
      printf("*** ERROR *** - No se reconoce la proyección elegida.");
      return -1;
    }
    if ( lon > pi ) {
      lon -= 2.0 * pi ;
    }
    printf("Vértice inferior derecho:\n" );
    fprintf( destinyFile, "%s", "Vértice inferior derecho:\n" );
    printf("  lon=%lf lat=%lf\n", lon*180.0/pi, lat*180/pi );
    fprintf( destinyFile, "  %s%lf %s%lf\n", "lon=", lon*180.0/pi, "lat=", lat*180/pi);

    /* Vértice superior izquierdo. */
    ix = 0;
    iy = 0;
    x = ( double ) ( ix - MAPA.X0 ) / MAPA.FACTORX + MAPA.Dx;
    y = ( double ) ( MAPA.Y0 - iy ) / MAPA.FACTORY + MAPA.Dy;
    // printf("x=%lf y=%lf\n", x, y);
    if ( MAPA.CPROY == 'O' )
      orthographic_inv( &lon, &lat, x, y, MAPA.XLON0, MAPA.XLAT0);
    else if ( MAPA.CPROY == 'M')
      mercator_inv( &lon, &lat, x, y, 0.0, 0.0 );
    else {
      printf("*** ERROR *** - No se reconoce la proyección elegida.");
      return -1;
    }    if ( lon > pi ) {
      lon -= 2.0 * pi;
    }
    printf("Vértice superior izquierdo:\n" );
    fprintf( destinyFile, "%s", "Vértice superior izquierdo:\n" );
    printf("  lon=%lf lat=%lf\n", lon*180.0/pi, lat*180/pi );
    fprintf( destinyFile, "  %s%lf %s%lf\n", "lon=", lon*180.0/pi, "lat=", lat*180/pi);

    /* Vértice superior derecho. */
    ix = MAPA.SIZEX - 1;
    iy = 0;
    x = ( double ) ( ix - MAPA.X0 ) / MAPA.FACTORX + MAPA.Dx;
    y = ( double ) ( MAPA.Y0 - iy ) / MAPA.FACTORY + MAPA.Dy;
    // printf("x=%lf y=%lf\n", x, y);
    if ( MAPA.CPROY == 'O' )
      orthographic_inv( &lon, &lat, x, y, MAPA.XLON0, MAPA.XLAT0);
    else if ( MAPA.CPROY == 'M')
      mercator_inv( &lon, &lat, x, y, 0.0, 0.0 );
    else {
      printf("*** ERROR *** - No se reconoce la proyección elegida.");
      return -1;
    }    if ( lon > pi ) {
      lon -= 2.0 * pi;
    }
    printf("Vértice superior derecho:\n" );
    fprintf( destinyFile, "%s", "Vértice superior derecho:\n" );
    printf("  lon=%lf lat=%lf\n", lon*180.0/pi, lat*180/pi );
    fprintf( destinyFile, "  %s%lf %s%lf\n", "lon=", lon*180.0/pi, "lat=", lat*180/pi);
    fprintf( destinyFile, "\n\n" );
  }

  fclose( destinyFile );
  return 0;
}