//#define TEST_MOLLWEIDE
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "sopnamsp.h"
#include "mollweide.h"

int main()
{
 FILE *file  = fopen("mollw.data","w");
 FILE *file2 = fopen("mollw2.data","w");
 int i,npoints=250000;
 double tvrai,lon,lat,x,y,tcalc,a;
 double longitude,latitude,tll,txy;

 tll=txy=0.;
 for(i=0;i<=npoints;i++) {
   //tvrai = 1.57078;
   tvrai = (double)i*(M_PI_2) /(double)npoints;
   a = 2.*tvrai+sin(2.*tvrai);
   lat = asin(a/M_PI);
   tcalc = MollWeide_XY(0.,lat,x,y);
   if(fabs(tcalc-tvrai)>tll) {tll = fabs(tcalc-tvrai); txy=tvrai;}
   if(file) fprintf(file,"%e %e %e %e\n",tvrai,tcalc,a,tcalc-tvrai);
 }
 printf("Maximum difference: %e for t=%.10e=Pi/2-%.10e\n",tll,txy,M_PI_2-txy);
 if(file) fclose(file);

 for(lon=0.;lon<360.;lon+=2.5) for(lat=-90.;lat<=90.1;lat+=2.5) {
   txy = MollWeide_XY(lon*M_PI/180.,lat*M_PI/180.,x,y);
   tll = MollWeide_LL(x,y,longitude,latitude);
   longitude *= 180./M_PI;
   latitude  *= 180./M_PI;
   if(file2) fprintf(file2,"%e %e %e %e %e %e %e %e %e %e\n"
     ,lon,lat,x,y,longitude,latitude,txy,tll,longitude-lon,latitude-lat);
 }
 if(file2) fclose(file2);
 return 0;
}