/* seviri_latlon.c */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include "geostat.h"

#define DEG2RAD M_PI/180.
#define RAD2DEG 180./M_PI

static float latitudes[3712*3712];
static float longitudes[3712*3712];

void seviri_compute_latlon(float *latitudes, float *longitudes) {
  const int nrows = 3712;
  const int ncols = 3712;
  const double NaN = 0./0.;
  int irow;
  int icol;
  double lat;
  double lon;
  geostat_err_t err;

  err = geostat_init("msg");
  assert(err == GEOSTAT_ERR_OK);
  
  
  for (irow = 0 ; irow < nrows ; irow++) {
    for (icol = 0 ; icol < ncols ; icol++) {
      
      err = geostat_lincol_to_latlon_double(irow + 1, icol + 1, &lat, &lon);
      if (err == GEOSTAT_ERR_OK) {
      	lat *= RAD2DEG;
      	lon *= RAD2DEG;
      }
	  else {
      	lat = NaN;
      	lon = NaN;
      }
      latitudes[irow*ncols+icol] = (float) lat;
      longitudes[irow*ncols+icol] = (float) lon;
    }
  }
}

int main(void) {

  const char file_latitudes[] = "seviri_latitudes.raw";
  const char file_longitudes[] = "seviri_longitudes.raw";
  FILE *stream;

  seviri_compute_latlon(latitudes, longitudes);

  stream = fopen(file_latitudes, "w");
  if (stream == NULL) {
    perror(file_latitudes);
    exit(EXIT_FAILURE);
  }

  fwrite(latitudes, 3712*3712*sizeof(float), 1, stream);

  if (fclose(stream) != 0) {
    perror(file_latitudes);
    exit(EXIT_FAILURE);
  }

  stream = fopen(file_longitudes, "w");
  if (stream == NULL) {
    perror(file_longitudes);
    exit(EXIT_FAILURE);
  }

  fwrite(longitudes, 3712*3712*sizeof(float), 1, stream);

  if (fclose(stream) != 0) {
    perror(file_longitudes);
    exit(EXIT_FAILURE);
  }
  
}