home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SRC.DI$ / YPRIME.C < prev    next >
Encoding:
C/C++ Source or Header  |  1993-03-26  |  2.0 KB  |  125 lines
  1. /*
  2.  
  3.   YPRIME.C    Sample .MEX file corresponding to YPRIME.M
  4.         Solves simple 3 body orbit problem 
  5.  
  6.   The calling syntax is:
  7.  
  8.             [yp] = yprime(t, y)
  9.  
  10.  
  11.   Marc Ullman   June 22, 1987
  12.   Copyright (c) 1987 The Mathworks, Inc.
  13.   All Righs Reserved
  14.  
  15.  
  16.   On PC platforms, must be compiled with /Alfw and /Gs option to generate
  17.   proper code.
  18.  
  19. */
  20.  
  21. #include <math.h>
  22. #include "mex.h"
  23.  
  24. /* Input Arguments */
  25.  
  26. #define    T_IN    prhs[0]
  27. #define    Y_IN    prhs[1]
  28.  
  29.  
  30. /* Output Arguments */
  31.  
  32. #define    YP_OUT    plhs[0]
  33.  
  34.  
  35. #define    max(A, B)    ((A) > (B) ? (A) : (B))
  36. #define    min(A, B)    ((A) < (B) ? (A) : (B))
  37.  
  38.  
  39. #define pi 3.14159265
  40.  
  41. static    double    mu = 1/82.45;
  42. static    double    mus = 1 - 1/82.45;
  43.  
  44.  
  45. static
  46. #ifdef __STDC__
  47. void yprime(
  48.     double    yp[],
  49.     double    *t,
  50.     double    y[]
  51.     )
  52. #else
  53. yprime(yp,t,y)
  54. double    yp[];
  55. double    *t,y[];
  56. #endif
  57. {
  58.     double    r1,r2;
  59.  
  60.     r1 = sqrt((y[0]+mu)*(y[0]+mu) + y[2]*y[2]);
  61.     r2 = sqrt((y[0]-mus)*(y[0]-mus) + y[2]*y[2]);
  62.  
  63.     yp[0] = y[1];
  64.     yp[1] = 2*y[3]+y[0]-mus*(y[0]+mu)/(r1*r1*r1)-mu*(y[0]-mus)/(r2*r2*r2);
  65.     yp[2] = y[3];
  66.     yp[3] = -2*y[1] + y[2] - mus*y[2]/(r1*r1*r1) - mu*y[2]/(r2*r2*r2);
  67.  
  68. }
  69.  
  70. #ifdef __STDC__
  71. void mexFunction(
  72.     int        nlhs,
  73.     Matrix    *plhs[],
  74.     int        nrhs,
  75.     Matrix    *prhs[]
  76.     )
  77. #else
  78. mexFunction(nlhs, plhs, nrhs, prhs)
  79. int nlhs, nrhs;
  80. Matrix *plhs[], *prhs[];
  81. #endif
  82. {
  83.     double    *yp;
  84.     double    *t,*y;
  85.     unsigned int    m,n;
  86.  
  87.     /* Check for proper number of arguments */
  88.  
  89.     if (nrhs != 2) {
  90.         mexErrMsgTxt("YPRIME requires two input arguments.");
  91.     } else if (nlhs > 1) {
  92.         mexErrMsgTxt("YPRIME requires one output argument.");
  93.     }
  94.  
  95.  
  96.     /* Check the dimensions of Y.  Y can be 4 X 1 or 1 X 4. */
  97.  
  98.     m = mxGetM(Y_IN);
  99.     n = mxGetN(Y_IN);
  100.     if (!mxIsNumeric(Y_IN) || mxIsComplex(Y_IN) || 
  101.         !mxIsFull(Y_IN)  || !mxIsDouble(Y_IN) ||
  102.         (max(m,n) != 4) || (min(m,n) != 1)) {
  103.         mexErrMsgTxt("YPRIME requires that Y be a 4 x 1 vector.");
  104.     }
  105.  
  106.  
  107.     /* Create a matrix for the return argument */
  108.  
  109.     YP_OUT = mxCreateFull(m, n, REAL);
  110.  
  111.  
  112.     /* Assign pointers to the various parameters */
  113.  
  114.     yp = mxGetPr(YP_OUT);
  115.  
  116.     t = mxGetPr(T_IN);
  117.     y = mxGetPr(Y_IN);
  118.  
  119.  
  120.     /* Do the actual computations in a subroutine */
  121.  
  122.     yprime(yp,t,y);
  123. }
  124.  
  125.