Want more features on Pastebin? Sign Up, it's FREE!
Guest

MATLAB MEX NRLMSISE BRIDGE

By: a guest on Nov 4th, 2011  |  syntax: C  |  size: 3.11 KB  |  views: 778  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. #include "mex.h"
  2. #include "nrlmsise-00.h"
  3. #include <math.h>
  4.  
  5. //returns atmospheric density, input is as follows:
  6. // [x, y, timePassed]
  7. //where x and y are points on a plane described below and
  8. //time is something I have yet to implement.
  9. void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] )
  10. {
  11.     //Declare the variables first because this is c
  12.   double x, y, time, theta, X, Y, Z, lat, lon, density;
  13.   struct nrlmsise_output output;
  14.   struct nrlmsise_input input;
  15.   struct nrlmsise_flags flags;
  16.   int startDay, i;
  17.   double *yp;
  18.   //Check to make sure the input is correct.
  19.  if(nrhs!=3)//Need three arguments
  20.   {
  21.     mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nrhs","Three inputs required.");
  22.   }
  23.   //returns one result, this appears to bu unnessicary.
  24.   /*if(nlhs!=1)Uncomment if you start experiencing errors.
  25.   {
  26.     mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nlhs","One output required.");
  27.   }*/
  28.   if(!mxIsDouble(prhs[0])||!mxIsDouble(prhs[1])||!mxIsDouble(prhs[2]))//They all need to be doubles
  29.   {
  30.     mexErrMsgIdAndTxt("MyToolbox:arrayProduct:notScalar","Wrong input");
  31.   }
  32.  
  33.   x = mxGetScalar(prhs[0]);//If you start using complex numbers, this might break
  34.   //double y;
  35.   y = mxGetScalar(prhs[1]);
  36.   time = mxGetScalar(prhs[2]);
  37.   //Input data is parsed, now convert x and y to lat, lon and alt:
  38.   //The x and y data will be converted to X, Y, Z, by saying it is on
  39.   //A plane,that crosses the XY plane at the X axis. x corisponds to X, while
  40.   //y coresponds to sqrt(Y^2+Z^2) and the angle between the xy plane and the XY
  41.   //plane will be represented as theta.  It is also convienient to say a y
  42.   //value of 0, which will correspond to a z value of 0, represents the
  43.   //intersection of the equator and the prime meridian.  This can be shifted
  44.   //by a time factor to increase accuracy if need be.
  45.   theta = .235;
  46.   startDay = 20; //The day of the year at which the satalite started falling
  47.   //Given x and y, z can be calculated as z=(Ax+By)/(-C)
  48.   //Converting sto circular coordinates
  49.   X=x;
  50.   Y=y*cos(theta);
  51.   Z=y*sin(theta);
  52.   lat = atan(Z/sqrt(X*X+Y*Y));
  53.   lon = atan2(Y,X);//Add in earth's rotation if you feel like it
  54.  
  55.    
  56.     i;
  57.     flags.switches[0]=0;
  58.     for(i=1;i<24;i++)
  59.         flags.switches[i]=1;
  60.     input.doy = 210;//Update for accuracy
  61.     input.year=0;//Not implemented
  62.     input.sec=29000;//Seconds passed, use input
  63.     input.alt = sqrt(x*x+y*y)-6380;//make sure Km are the correct units
  64.     input.g_lat=lat*180/3.14159;//latitude in degrees
  65.     input.g_long=lon*180/3.14159;//longitude in degrees
  66.     input.lst = input.sec/3600+input.g_long/15;//random junk
  67.     input.f107A=150;//not important
  68.     input.f107=150;//doesn't matter
  69.     input.ap=4;//has little effect
  70.     gtd7(&input, &flags, &output);
  71.     //for debuging, uncomment following line
  72.    // 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]);
  73.     plhs[0] = mxCreateDoubleScalar(output.d[5]);
  74. }
  75.  
  76.  
clone this paste RAW Paste Data