/*=================================================================
 *
 * YPRIME.C	Sample .MEX file corresponding to YPRIME.M
 *	        Solves simple 3 body orbit problem 
 *
 * The calling syntax is:
 *
 *		[yp] = yprime(t, y)
 *
 *  You may also want to look at the corresponding M-code, yprime.m.
 *
 * This is a MEX-file for MATLAB.  
 * Copyright 1984-2006 The MathWorks, Inc.
 *
 *=================================================================*/
/* $Revision: 1.10.6.4 $ */
#include <math.h>
#include "mex.h"


static void GEOM_OUT_C( 
        double Ao[],
        double Bo[],
        double Lo[],
        double Vo[],
        double Yo[],
        double Ax[],
        double Ay[],
        double Bx[],
        double By[])
{   
    int i;
    double theta, pi = 3.14159265359;
    double tmp_dax, tmp_day;
    double tmp_dbx, tmp_dby;
    double tmp_dlx, tmp_dly;
    
    double tmp;
    
    for (i=0;i<80;i++){
        Vo[i] = 0;
    }
    Yo[0] = 0;
    
    for (i=0;i<80;i++){
        /* Apical surface */
        tmp_dax = Ax[(i+1)%80] - Ax[i];
        tmp_day = Ay[(i+1)%80] - Ay[i];
        Ao[i] = sqrt(tmp_dax*tmp_dax + tmp_day*tmp_day);
               
        /* Basal surface */
        tmp_dbx = Bx[(i+1)%80] - Bx[i];
        tmp_dby = By[(i+1)%80] - By[i];
        Bo[i] = sqrt(tmp_dbx*tmp_dbx + tmp_dby*tmp_dby);
        
        /* Lateral surface */
        tmp_dlx = Bx[(i+1)%80] - Ax[(i+1)%80];
        tmp_dly = By[(i+1)%80] - Ay[(i+1)%80];
        Lo[i] = sqrt(tmp_dlx*tmp_dlx + tmp_dly*tmp_dly);
        
        /* Sum up the volume */
        
        tmp = (Ax[(i+1)%80] + Bx[(i+1)%80])*tmp_dly/4.0 - (Ay[(i+1)%80] + By[(i+1)%80])*tmp_dlx/4.0;
        Vo[(i+1)%80] -= tmp;
        Vo[i] += tmp;
        Vo[i] += (Ax[(i+1)%80] + Ax[i])*tmp_day/4.0 - (Ay[(i+1)%80] + Ay[i])*tmp_dax/4.0;
        Vo[i] -= (Bx[(i+1)%80] + Bx[i])*tmp_dby/4.0 - (By[(i+1)%80] + By[i])*tmp_dbx/4.0;
        
        /* Sum up the yolk */
        Yo[0] += (Bx[(i+1)%80] + Bx[i])*tmp_dby/4.0 - (By[(i+1)%80] + By[i])*tmp_dbx/4.0;
   }
    
    return;
}

void mexFunction( int nlhs, mxArray *plhs[], 
		  int nrhs, const mxArray *prhs[] )
     
{ 
    double *Ax,*Ay,*Bx,*By;
    double *Ao,*Bo,*Lo,*Vo,*Yo;
    mwSize m,n; 
    
    /* Check for proper number of arguments */
    
    if (nrhs != 4) { 
	mexErrMsgTxt("Incorrect number of input arguments."); 
    } 
    if (nlhs != 5) {
	mexErrMsgTxt("Incorrect number of output arguments."); 
    } 
    
    /* Create a matrix for the return argument */ 
    plhs[0] = mxCreateDoubleMatrix(1, 80, mxREAL);
    plhs[1] = mxCreateDoubleMatrix(1, 80, mxREAL);
    plhs[2] = mxCreateDoubleMatrix(1, 80, mxREAL);
    plhs[3] = mxCreateDoubleMatrix(1, 80, mxREAL);
    plhs[4] = mxCreateDoubleMatrix(1, 1, mxREAL);
    
    /* Assign pointers to the various parameters */ 

    Ao = mxGetPr(plhs[0]);
    Bo = mxGetPr(plhs[1]);
    Lo = mxGetPr(plhs[2]);
    Vo = mxGetPr(plhs[3]);
    Yo = mxGetPr(plhs[4]);

    Ax = mxGetPr(prhs[0]); 
    Ay = mxGetPr(prhs[1]);
    Bx = mxGetPr(prhs[2]); 
    By = mxGetPr(prhs[3]);
    
    /* Do the actual computations in a subroutine */
    GEOM_OUT_C(Ao,Bo,Lo,Vo,Yo,Ax,Ay,Bx,By); 
    return;
    
}


