rhotheta2xy_gb

C/C++ source

0001 /*
0002  rhotheta2xy_gb.c : MEX file giving out the corresponding cartesian couple from an input curvilinear couple
0003  
0004  INPUT:
0005  - Magnetic equilibrium                           : equil_IN
0006  - Ray parameters                                 : rayparam_IN
0007  - Equilibrium parameters                         : equilparam_IN
0008  - Radial coordinate                              : rho_IN
0009  - Poloidal coordinate                            : theta_IN
0010  
0011  OUTPUT:
0012  - Cartesian coordinate x : x_OUT
0013  - Cartesian coordinate y : y_OUT
0014  
0015  By Guillaume Brochard (CEA-DRFC, guillaume.brochard@cea.fr)
0016 */
0017 #define DUMP(u) mexPrintf("%s = %e\n", #u, u)
0018 
0019 #include <math.h>
0020 #include "mex.h"
0021 #include <complex.h>
0022 #include "structures_raytracing_yp.h"
0023 #include <time.h>
0024 
0025 /*input arguments*/
0026 
0027 #define    equil_IN      prhs[0]
0028 #define rayparam_IN   prhs[1]
0029 #define equilparam_IN prhs[2]
0030 #define rho_IN       prhs[3]
0031 #define theta_IN      prhs[4]
0032 
0033 /*output arguments*/
0034 
0035 #define x_OUT       plhs[0]
0036 #define y_OUT     plhs[1]
0037 
0038 void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
0039 {
0040     double *x,*y,*rho,*theta;
0041     
0042     rho = mxCalloc(1,sizeof(double));
0043     theta = mxCalloc(1,sizeof(double));
0044     
0045     /*Declaration of the input parameters*/
0046     
0047     rho         =   mxGetPr(rho_IN);
0048     theta       =   mxGetPr(theta_IN);
0049 
0050     
0051     /*Declaration of the output parameters*/
0052     
0053     x_OUT = mxCreateDoubleScalar(mxREAL);
0054     y_OUT = mxCreateDoubleScalar(mxREAL);
0055     x = mxGetPr(x_OUT);
0056     y = mxGetPr(y_OUT);
0057     
0058     
0059     // Declaration of the structures
0060     
0061     rayparam_format rayparam;
0062     rayequil_format rayequil;
0063     equilparam_format equilparam;
0064     
0065     // Loading of the magnetic equilibrium
0066     
0067     rayequil.x_fit       = load_grid_2D(equil_IN,"x_fit",0);
0068     rayequil.y_fit       = load_grid_2D(equil_IN,"y_fit",0);
0069     rayparam             = load_ray_param(rayparam_IN);
0070     equilparam           = load_equil_param(equilparam_IN);
0071     
0072     //Scheme
0073     
0074     double a,b,dummy;
0075     ppvald_2D(rayequil.x_fit,rayparam.dS[0],rho[0],theta[0],x,&dummy,&dummy,&dummy,&dummy);//Gives out x
0076     ppvald_2D(rayequil.y_fit,rayparam.dS[0],rho[0],theta[0],y,&dummy,&dummy,&dummy,&dummy);//Gives out y
0077     x[0]=x[0]/equilparam.ap[0];
0078     y[0]=y[0]/equilparam.ap[0];
0079     return;
0080 }