Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #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]);
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement