#include "mex.h"
#include "nrlmsise-00.h"
#include <math.h>

//returns atmospheric density, input is as follows:
// [x, y, timePassed]
//where x and y are points on a plane described below and 
//time is something I have yet to implement.
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] )
{
    //Declare the variables first because this is c
  double x, y, time, theta, X, Y, Z, lat, lon, density;
  struct nrlmsise_output output;
  struct nrlmsise_input input;
  struct nrlmsise_flags flags;
  int startDay, i;
  double *yp;
  //Check to make sure the input is correct.
 if(nrhs!=3)//Need three arguments
  {
    mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nrhs","Three inputs required.");
  }
  //returns one result, this appears to bu unnessicary.
  /*if(nlhs!=1)Uncomment if you start experiencing errors.
  {
    mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nlhs","One output required.");
  }*/
  if(!mxIsDouble(prhs[0])||!mxIsDouble(prhs[1])||!mxIsDouble(prhs[2]))//They all need to be doubles
  {
    mexErrMsgIdAndTxt("MyToolbox:arrayProduct:notScalar","Wrong input");
  }
  
  x = mxGetScalar(prhs[0]);//If you start using complex numbers, this might break
  //double y;
  y = mxGetScalar(prhs[1]);
  time = mxGetScalar(prhs[2]);
  //Input data is parsed, now convert x and y to lat, lon and alt:
  //The x and y data will be converted to X, Y, Z, by saying it is on
  //A plane,that crosses the XY plane at the X axis. x corisponds to X, while
  //y coresponds to sqrt(Y^2+Z^2) and the angle between the xy plane and the XY
  //plane will be represented as theta.  It is also convienient to say a y 
  //value of 0, which will correspond to a z value of 0, represents the 
  //intersection of the equator and the prime meridian.  This can be shifted
  //by a time factor to increase accuracy if need be.
  theta = .235;
  startDay = 20; //The day of the year at which the satalite started falling
  //Given x and y, z can be calculated as z=(Ax+By)/(-C)
  //Converting sto circular coordinates
  X=x;
  Y=y*cos(theta);
  Z=y*sin(theta);
  lat = atan(Z/sqrt(X*X+Y*Y));
  lon = atan2(Y,X);//Add in earth's rotation if you feel like it
  
    
    i;
    flags.switches[0]=0;
    for(i=1;i<24;i++)
        flags.switches[i]=1;
    input.doy = 210;//Update for accuracy
    input.year=0;//Not implemented
    input.sec=29000;//Seconds passed, use input
    input.alt = sqrt(x*x+y*y)-6380;//make sure Km are the correct units
    input.g_lat=lat*180/3.14159;//latitude in degrees
    input.g_long=lon*180/3.14159;//longitude in degrees
    input.lst = input.sec/3600+input.g_long/15;//random junk
    input.f107A=150;//not important
    input.f107=150;//doesn't matter
    input.ap=4;//has little effect
    gtd7(&input, &flags, &output);
    //for debuging, uncomment following line
   // mexPrintf("\nCalling function with following parameters:\n\tDay of year: %f\n\tSeconds: %f\n\tAltitude: %f\n\tLatitude: %f\n\tLongitude: %f\nThe density is %Ef g/cm^3\n",input.doy,input.sec,input.alt, input.g_lat, input.g_long,output.d[5]);
    plhs[0] = mxCreateDoubleScalar(output.d[5]);
}