/*=================================================================
 *
 * 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 INIT_C( 
        double Ax[],
        double Ay[],
        double Bx[],
        double By[],
        double Ao_con[],
        double Ao[],
        double Ka_lin[],
        double Ka_con[],
        double Ka[],
        double Bo[],
        double Kb_lin[],
        double Kb[],
        double Lo[],
        double Kl[],
        double Vo[],
        double Kv[],
        double Yo[],
        double Ky[],
        double R_out[],
        double R_in[],
        double alpha[],
        double S[],
        double KA_val[],
        double KB_val[],
        double KL_val[],
        double KV_val[],
        double KY_val[])
{   
    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++){
        theta = -pi/2.0 + i*2.0*pi/80.0;
        Ax[i] = R_out[0]*cos(theta);
        Ay[i] = R_out[0]*sin(theta);
        Bx[i] = R_in[0]*cos(theta);
        By[i] = R_in[0]*sin(theta);
        
        theta = (-pi*R_out[0] + i*2*pi*R_out[0]/79.0)/S[0];
        Ka_lin[i] = alpha[0]*exp(-theta*theta/2);
        
        theta = (-pi*R_out[0] + i*2*pi*R_out[0]/79.0)/S[1];
        Ka_con[i] = alpha[1]*exp(-theta*theta/2);
        
        if (S[2] > 0){
        	theta = (-pi*R_out[0] + i*2*pi*R_out[0]/79.0)/S[2];
        	Kb_lin[i] = alpha[2]*exp(-theta*theta/2);
        }
        else {
        	Kb_lin[i] = alpha[2];
        }
        Ka[i] = KA_val[0];
        Kl[i] = KL_val[0];
        Kb[i] = KB_val[0];
        Kv[i] = KV_val[0];
        
        Vo[i] = 0;
    }
    Ky[0] = KY_val[0];Ky[1] = KY_val[1];
    Yo[0] = 0;
    
    for (i=0;i<80;i++){
        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);
        Ao_con[i] = sqrt(tmp_dax*tmp_dax + tmp_day*tmp_day);
               
        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);
        
        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);
        
        
        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;
        
        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 *R_out,*R_in,*alpha,*S,*KA_val,*KB_val,*KL_val,*KV_val,*KY_val; 
    double *Ax,*Ay,*Bx,*By,*Ao_con,*Ao,*Ka_lin,*Ka_con,*Ka,*Bo,*Kb_lin,*Kb,*Lo,*Kl,*Vo,*Kv,*Yo,*Ky;
    mwSize m,n; 
    
    /* Check for proper number of arguments */
    
    if (nrhs != 9) { 
	mexErrMsgTxt("Incorrect number of input arguments."); 
    } 
    if (nlhs != 18) {
	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, 80, mxREAL);
    plhs[5] = mxCreateDoubleMatrix(1, 80, mxREAL);
    plhs[6] = mxCreateDoubleMatrix(1, 80, mxREAL);
    plhs[7] = mxCreateDoubleMatrix(1, 80, mxREAL);
    plhs[8] = mxCreateDoubleMatrix(1, 80, mxREAL);
    plhs[9] = mxCreateDoubleMatrix(1, 80, mxREAL);
    plhs[10] = mxCreateDoubleMatrix(1, 80, mxREAL);
    plhs[11] = mxCreateDoubleMatrix(1, 80, mxREAL);
    plhs[12] = mxCreateDoubleMatrix(1, 80, mxREAL);
    plhs[13] = mxCreateDoubleMatrix(1, 80, mxREAL);
    plhs[14] = mxCreateDoubleMatrix(1, 80, mxREAL);
    plhs[15] = mxCreateDoubleMatrix(1, 80, mxREAL);
    plhs[16] = mxCreateDoubleMatrix(1, 1, mxREAL);
    plhs[17] = mxCreateDoubleMatrix(1, 2, mxREAL);
    
    
    /* Assign pointers to the various parameters */ 
    Ax = mxGetPr(plhs[0]);
    Ay = mxGetPr(plhs[1]);
    Bx = mxGetPr(plhs[2]);
    By = mxGetPr(plhs[3]);
    Ao_con = mxGetPr(plhs[4]);
    Ao = mxGetPr(plhs[5]);
    Ka_lin = mxGetPr(plhs[6]);
    Ka_con = mxGetPr(plhs[7]);
    Ka = mxGetPr(plhs[8]);
    Bo = mxGetPr(plhs[9]);
    Kb_lin = mxGetPr(plhs[10]);
    Kb = mxGetPr(plhs[11]);
    Lo = mxGetPr(plhs[12]);
    Kl = mxGetPr(plhs[13]);
    Vo = mxGetPr(plhs[14]);
    Kv = mxGetPr(plhs[15]);
    Yo = mxGetPr(plhs[16]);
    Ky = mxGetPr(plhs[17]);
    
    R_out = mxGetPr(prhs[0]); 
    R_in = mxGetPr(prhs[1]);
    alpha = mxGetPr(prhs[2]); 
    S = mxGetPr(prhs[3]);
    KA_val = mxGetPr(prhs[4]);
    KB_val = mxGetPr(prhs[5]);
    KL_val = mxGetPr(prhs[6]);
    KV_val = mxGetPr(prhs[7]);
    KY_val = mxGetPr(prhs[8]);
    
    /* Do the actual computations in a subroutine */
    INIT_C(Ax,Ay,Bx,By,Ao_con,Ao,Ka_lin,Ka_con,Ka,Bo,Kb_lin,Kb,Lo,Kl,Vo,Kv,Yo,Ky,R_out,R_in,alpha,S,KA_val,KB_val,KL_val,KV_val,KY_val); 
    return;
    
}


