#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]);
}