test_grid_substitution

C/C++ source

0001 /*
0002     test_grid_substitution_gb.c    MEX file the mean precision, execution time and dispersion of the (x,y)->(rho,theta)->(x,y) on a chosen grid
0003  
0004     test_grid_substitution_gb(equil_IN,rayparam_IN,equilparam_IN,ngrid_IN,length_IN)
0005  
0006     INPUT:
0007  - Magnetic equilibrium                           : equil_IN
0008  - Ray parameters                                 : rayparam_IN
0009  - Equilibrium parameters                         : equilparam_IN
0010  - Number of points on a side of the grid         : ngrid_IN
0011  - Length of the grid                             : length_IN
0012  
0013     OUTPUT:
0014     
0015  - Three dimension array : data
0016  - Mean precision        : data[0]
0017  - Mean dispersion       : data[1]
0018  - Mean execution time   : data[2]
0019  
0020  By Guillaume Brochard (CEA-DRFC, guillaume.brochard@cea.fr)
0021  */
0022 
0023 #define DUMP(u) mexPrintf("%s = %e\n", #u, u)
0024 
0025 #include <math.h>
0026 #include "mex.h"
0027 #include <complex.h>
0028 #include "structures_raytracing_yp.h"
0029 #include <time.h>
0030 
0031 /*input arguments*/
0032 
0033 #define    equil_IN        prhs[0]
0034 #define rayparam_IN     prhs[1]
0035 #define equilparam_IN   prhs[2]
0036 #define ngrid_IN        prhs[3]
0037 #define length_IN       prhs[4]
0038 
0039 #define data_OUT       plhs[0]
0040 
0041 void rho1D();
0042 void rho1D_newton();
0043 void substitution();
0044 void substitution_newton();
0045 
0046 void rho1D(rayequil_format rayequil,rayparam_format rayparam, equilparam_format equilparam, double *rho,double *theta,double *Y)
0047 {
0048     double down[2],up[2],res[2];
0049     down[0]=0;
0050     up[0]=1;
0051     double dummy,x,y,dxdrho,dydrho;
0052     int k = 0;
0053     res[1]=1;
0054     int kmax = 60;
0055     double prec = 1e-10;
0056     while(k < kmax && fabs(res[1]) > prec)
0057     {
0058         k = k +1;
0059         
0060         res[0] = (down[0] + up[0])/2;
0061         
0062         ppvald_2D(rayequil.x_fit,rayparam.dS[0],down[0],theta[0],&x,&dxdrho,&dummy,&dummy,&dummy,&dummy,&dummy);
0063         ppvald_2D(rayequil.y_fit,rayparam.dS[0],down[0],theta[0],&y,&dydrho,&dummy,&dummy,&dummy,&dummy,&dummy);
0064         x=x/equilparam.ap[0];//Mandatory normalization
0065         y=y/equilparam.ap[0];//Mandatory normalization
0066         dxdrho=dxdrho/equilparam.ap[0];//Mandatory normalization
0067         dydrho=dydrho/equilparam.ap[0];//Mandatory normalization
0068         down[1] = dxdrho*(x-Y[0]) + dydrho*(y - Y[1]);
0069         
0070         ppvald_2D(rayequil.x_fit,rayparam.dS[0],res[0],theta[0],&x,&dxdrho,&dummy,&dummy,&dummy,&dummy,&dummy);
0071         ppvald_2D(rayequil.y_fit,rayparam.dS[0],res[0],theta[0],&y,&dydrho,&dummy,&dummy,&dummy,&dummy,&dummy);
0072         x=x/equilparam.ap[0];//Mandatory normalization
0073         y=y/equilparam.ap[0];//Mandatory normalization
0074         dxdrho=dxdrho/equilparam.ap[0];//Mandatory normalization
0075         dydrho=dydrho/equilparam.ap[0];//Mandatory normalization
0076         res[1] = dxdrho*(x-Y[0]) + dydrho*(y - Y[1]);
0077         
0078         if(down[1]*res[1]>0)
0079         {
0080             down[0]=res[0];
0081         }
0082         else
0083         {
0084             up[0]=res[0];
0085         }
0086     }
0087     rho[0]=res[0];
0088 }
0089 
0090 void substitution(rayequil_format rayequil,rayparam_format rayparam, equilparam_format equilparam, double *rho,double *theta,double *Y)
0091 {
0092     theta[0] = atan2(Y[1],Y[0]);
0093     double down_t[2],up_t[2],res_t[2];
0094     res_t[0] = theta[0];
0095     down_t[0] = theta[0]*0.9;
0096     up_t[0] = theta[0]*1.1;
0097     double dummy,x,y,dxdtheta,dydtheta,rho_d,rho_m;
0098     int k=0;
0099     int kmax = 60;
0100     double prec = 1e-11;
0101     res_t[1] = 1;
0102     while(k < kmax && fabs(res_t[1]) > prec)
0103     {
0104         k = k + 1;
0105         res_t[0] = (down_t[0] + up_t[0])/2;
0106         rho1D(rayequil,rayparam,equilparam,&rho_d,&down_t[0],Y);
0107         rho1D(rayequil,rayparam,equilparam,&rho_m,&res_t[0],Y);
0108         
0109         ppvald_2D(rayequil.x_fit,rayparam.dS[0],rho_d,down_t[0],&x,&dummy,&dxdtheta,&dummy,&dummy,&dummy,&dummy);
0110         ppvald_2D(rayequil.y_fit,rayparam.dS[0],rho_d,down_t[0],&y,&dummy,&dydtheta,&dummy,&dummy,&dummy,&dummy);
0111         x=x/equilparam.ap[0];//Mandatory normalization
0112         y=y/equilparam.ap[0];//Mandatory normalization
0113         dxdtheta=dxdtheta/equilparam.ap[0];//Mandatory normalization
0114         dydtheta=dydtheta/equilparam.ap[0];//Mandatory normalization
0115         down_t[1] = dxdtheta*(x-Y[0]) + dydtheta*(y-Y[1]);
0116         
0117         ppvald_2D(rayequil.x_fit,rayparam.dS[0],rho_m,res_t[0],&x,&dummy,&dxdtheta,&dummy,&dummy,&dummy,&dummy);
0118         ppvald_2D(rayequil.y_fit,rayparam.dS[0],rho_m,res_t[0],&y,&dummy,&dydtheta,&dummy,&dummy,&dummy,&dummy);
0119         x=x/equilparam.ap[0];//Mandatory normalization
0120         y=y/equilparam.ap[0];//Mandatory normalization
0121         dxdtheta=dxdtheta/equilparam.ap[0];//Mandatory normalization
0122         dydtheta=dydtheta/equilparam.ap[0];//Mandatory normalization
0123         res_t[1] = dxdtheta*(x-Y[0]) + dydtheta*(y-Y[1]);
0124         
0125         if(down_t[1]*res_t[1]>0)
0126         {
0127             down_t[0] = res_t[0];
0128         }
0129         else
0130         {
0131             up_t[0] = res_t[0];
0132         }
0133     }
0134     theta[0] = res_t[0];
0135     rho1D(rayequil,rayparam,equilparam,rho,theta,Y);
0136 }
0137 
0138 void rho1D_newton(rayequil_format rayequil,rayparam_format rayparam, equilparam_format equilparam, double *rho,double *theta,double *Y)
0139 {
0140     double down[2],res[2];
0141     down[0]=0.01;
0142     res[0]=0.01000000000001;
0143     double dummy,x,y,dxdrho,dydrho;
0144     int k = 0;
0145     res[1]=1;
0146     int kmax = 10;
0147     double prec = 1e-10;
0148     double save1,save2;
0149     while(k < kmax && fabs(res[1]) > prec)
0150     {
0151         ppvald_2D(rayequil.x_fit,rayparam.dS[0],down[0],theta[0],&x,&dxdrho,&dummy,&dummy,&dummy,&dummy,&dummy);
0152         ppvald_2D(rayequil.y_fit,rayparam.dS[0],down[0],theta[0],&y,&dydrho,&dummy,&dummy,&dummy,&dummy,&dummy);
0153         
0154         x=x/equilparam.ap[0];//Mandatory normalization
0155         y=y/equilparam.ap[0];//Mandatory normalization
0156         dxdrho=dxdrho/equilparam.ap[0];//Mandatory normalization
0157         dydrho=dydrho/equilparam.ap[0];//Mandatory normalization
0158         
0159         down[1] = dxdrho*(x-Y[0]) + dydrho*(y - Y[1]);
0160         
0161         ppvald_2D(rayequil.x_fit,rayparam.dS[0],res[0],theta[0],&x,&dxdrho,&dummy,&dummy,&dummy,&dummy,&dummy);
0162         ppvald_2D(rayequil.y_fit,rayparam.dS[0],res[0],theta[0],&y,&dydrho,&dummy,&dummy,&dummy,&dummy,&dummy);
0163         
0164         x=x/equilparam.ap[0];//Mandatory normalization
0165         y=y/equilparam.ap[0];//Mandatory normalization
0166         dxdrho=dxdrho/equilparam.ap[0];//Mandatory normalization
0167         dydrho=dydrho/equilparam.ap[0];//Mandatory normalization
0168         
0169         res[1] = dxdrho*(x-Y[0]) + dydrho*(y - Y[1]);
0170         
0171         save1 = res[0];
0172         save2 = res[1];
0173         res[0] = res[0] - res[1]*((res[0] - down[0])/(res[1] - down[1]));
0174         down[0] = save1;
0175         down[1] = save2;
0176         
0177         k = k +1;
0178         
0179     }
0180     rho[0]=res[0];
0181 }
0182 
0183 void substitution_newton(rayequil_format rayequil,rayparam_format rayparam, equilparam_format equilparam, double *rho,double *theta,double *Y)
0184 {
0185     theta[0] = atan2(Y[1],Y[0]);
0186     double down_t[2],res_t[2];
0187     double dummy,x,y,dxdtheta,dydtheta,rho_d,rho_m;
0188     down_t[0] = theta[0];
0189     
0190     rho1D_newton(rayequil,rayparam,equilparam,&rho_d,&down_t[0],Y);
0191     ppvald_2D(rayequil.x_fit,rayparam.dS[0],rho_d,down_t[0],&x,&dummy,&dxdtheta,&dummy,&dummy,&dummy,&dummy);
0192     ppvald_2D(rayequil.y_fit,rayparam.dS[0],rho_d,down_t[0],&y,&dummy,&dydtheta,&dummy,&dummy,&dummy,&dummy);
0193     x=x/equilparam.ap[0];//Mandatory normalization
0194     y=y/equilparam.ap[0];//Mandatory normalization
0195     dxdtheta=dxdtheta/equilparam.ap[0];//Mandatory normalization
0196     dydtheta=dydtheta/equilparam.ap[0];//Mandatory normalization
0197     down_t[1] = dxdtheta*(x-Y[0]) + dydtheta*(y-Y[1]);
0198     
0199     res_t[0] = 1.0000000000001*down_t[0];
0200     int kt=0;
0201     int ktmax = 10;
0202     double prec = 1e-11;
0203     res_t[1] = 1;
0204     double save1,save2;
0205     while(kt < ktmax && fabs(res_t[1]) > prec)
0206     {
0207         rho1D_newton(rayequil,rayparam,equilparam,&rho_m,&res_t[0],Y);
0208         ppvald_2D(rayequil.x_fit,rayparam.dS[0],rho_m,res_t[0],&x,&dummy,&dxdtheta,&dummy,&dummy,&dummy,&dummy);
0209         ppvald_2D(rayequil.y_fit,rayparam.dS[0],rho_m,res_t[0],&y,&dummy,&dydtheta,&dummy,&dummy,&dummy,&dummy);
0210         x=x/equilparam.ap[0];//Mandatory normalization
0211         y=y/equilparam.ap[0];//Mandatory normalization
0212         dxdtheta=dxdtheta/equilparam.ap[0];//Mandatory normalization
0213         dydtheta=dydtheta/equilparam.ap[0];//Mandatory normalization
0214         res_t[1] = dxdtheta*(x-Y[0]) + dydtheta*(y-Y[1]);
0215         
0216         save1 = res_t[0];
0217         save2 = res_t[1];
0218         res_t[0] = res_t[0] - res_t[1]*((res_t[0] - down_t[0])/(res_t[1] - down_t[1]));
0219         down_t[0] = save1;
0220         down_t[1] = save2;
0221         
0222         kt = kt + 1;
0223         
0224     }
0225     theta[0] = res_t[0];
0226     rho1D_newton(rayequil,rayparam,equilparam,rho,theta,Y);
0227 }
0228 
0229 void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
0230 {
0231     double *data,*Y,*rho,*theta,*ngrid,*length;
0232     
0233     rho = mxCalloc(1,sizeof(double));
0234     theta = mxCalloc(1,sizeof(double));
0235     Y = mxCalloc(4,sizeof(double));
0236     ngrid = mxCalloc(1,sizeof(double));
0237     length = mxCalloc(1,sizeof(double));
0238     
0239     
0240     /*Declaration of the input parameters*/
0241     
0242     ngrid          =   mxGetPr(ngrid_IN);
0243     length         =   mxGetPr(length_IN);
0244     
0245     // Check for proper number of arguments
0246     
0247     if (!mxIsStruct(prhs[0])) mexErrMsgTxt("1st input argument of call.mex must be a structure (equil_fit)");
0248     if (nrhs != 5) mexErrMsgTxt("Wrong number of input arguments for raytracing_yp.mex");// The number of input arguments is fixed
0249 
0250     // Declaration of the structures
0251     
0252     rayparam_format rayparam;
0253     rayequil_format rayequil;
0254     equilparam_format equilparam;
0255 
0256     // Loading of the magnetic equilibrium
0257     
0258     rayequil.x_fit       = load_grid_2D(equil_IN,"x_fit",0);
0259     rayequil.y_fit       = load_grid_2D(equil_IN,"y_fit",0);
0260     rayparam             = load_ray_param(rayparam_IN);
0261     equilparam           = load_equil_param(equilparam_IN);
0262 
0263     // Scheme
0264     
0265     double dummy;
0266     int n = (int) ngrid[0];//Number of grid points
0267     data_OUT = mxCreateDoubleMatrix(1,3,mxREAL);
0268     data = mxGetPr(data_OUT);
0269     double *x0;
0270     double *y0;
0271     x0 = mxCalloc(1,sizeof(double));
0272     y0 = mxCalloc(1,sizeof(double));
0273     x0[0] = -length[0]/2;
0274     y0[0] = -length[0]/2;
0275     double N[n][n];
0276     double meand = 0;
0277     double dispd = 0;
0278     double dtd = 0;
0279     clock_t start,end;
0280     double cpu;
0281     int i,j,l;
0282     double a,b;
0283     meand = 0;
0284     dispd = 0;
0285     dtd = 0;
0286 
0287     for(i=0;ifor(j=0;jreturn;
0329 }