xy2rhotheta_gb

PURPOSE ^

C/C++ source

SYNOPSIS ^

C/C++ source

DESCRIPTION ^

C/C++ source

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 /*
0002  rhotheta2xy_gb.c : MEX file giving out the corresponding curvilinear couple from an input cartesian couple
0003  
0004  INPUT:
0005  - Magnetic equilibrium                           : equil_IN
0006  - Ray parameters                                 : rayparam_IN
0007  - Equilibrium parameters                         : equilparam_IN
0008  - Cartesian coordinate x                         : x_IN
0009  - Cartesian coordinate y                         : y_IN
0010  
0011  OUTPUT:
0012  - Radial coordinate   : rho_OUT
0013  - Poloidal coordinate : theta_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 x_IN       prhs[3]
0031 #define y_IN      prhs[4]
0032 
0033 /*output arguments*/
0034 
0035 #define rho_OUT       plhs[0]
0036 #define theta_OUT     plhs[1]
0037 
0038 void rho1D_dichotomy();
0039 void rho1D_newton_numeric();
0040 void rho1D_newton_analytic();
0041 void substitution_dichotomy_2D();
0042 void substitution_newton_numeric_2D();
0043 void substitution_newton_analytic_2D();
0044 
0045 void rho1D_dichotomy(rayequil_format rayequil,rayparam_format rayparam, equilparam_format equilparam, double *rho,double *theta,double *Y)
0046 {
0047     double down[2],up[2],res[2];
0048     down[0]=0;
0049     up[0]=1;
0050     double dummy,x,y,dxdrho,dydrho;
0051     int k = 0;
0052     res[1]=1;
0053     int kmax = 60;
0054     double prec = 1e-10;
0055     while(k < kmax && fabs(res[1]) > prec)
0056     {
0057         k = k +1;
0058         
0059         res[0] = (down[0] + up[0])/2;
0060         
0061         ppvald_2D(rayequil.x_fit,rayparam.dS[0],down[0],theta[0],&x,&dxdrho,&dummy,&dummy,&dummy,&dummy);
0062         ppvald_2D(rayequil.y_fit,rayparam.dS[0],down[0],theta[0],&y,&dydrho,&dummy,&dummy,&dummy,&dummy);
0063         x=x/equilparam.ap[0];//Mandatory normalization
0064         y=y/equilparam.ap[0];//Mandatory normalization
0065         dxdrho=dxdrho/equilparam.ap[0];//Mandatory normalization
0066         dydrho=dydrho/equilparam.ap[0];//Mandatory normalization
0067         down[1] = dxdrho*(x-Y[0]) + dydrho*(y - Y[1]);
0068         
0069         ppvald_2D(rayequil.x_fit,rayparam.dS[0],res[0],theta[0],&x,&dxdrho,&dummy,&dummy,&dummy,&dummy);
0070         ppvald_2D(rayequil.y_fit,rayparam.dS[0],res[0],theta[0],&y,&dydrho,&dummy,&dummy,&dummy,&dummy);
0071         x=x/equilparam.ap[0];//Mandatory normalization
0072         y=y/equilparam.ap[0];//Mandatory normalization
0073         dxdrho=dxdrho/equilparam.ap[0];//Mandatory normalization
0074         dydrho=dydrho/equilparam.ap[0];//Mandatory normalization
0075         res[1] = dxdrho*(x-Y[0]) + dydrho*(y - Y[1]);
0076         
0077         if(down[1]*res[1]>0)
0078         {
0079             down[0]=res[0];
0080         }
0081         else
0082         {
0083             up[0]=res[0];
0084         }
0085     }
0086     rho[0]=res[0];
0087 }
0088 
0089 void rho1D_newton_numeric(rayequil_format rayequil,rayparam_format rayparam, equilparam_format equilparam, double *rho,double *theta,double *Y)
0090 {
0091     double down[2],res[2];
0092     down[0]=0.05;
0093     res[0]=0.050000001;
0094     double dummy,x,y,dxdrho,dydrho;
0095     int k = 0;
0096     res[1]=1;
0097     int kmax = 10;
0098     double prec = 1e-9;
0099     double save;
0100     while(k < kmax && fabs(res[1]) > prec)
0101     {
0102         
0103         ppvald_2D(rayequil.x_fit,rayparam.dS[0],down[0],theta[0],&x,&dxdrho,&dummy,&dummy,&dummy,&dummy);
0104         ppvald_2D(rayequil.y_fit,rayparam.dS[0],down[0],theta[0],&y,&dydrho,&dummy,&dummy,&dummy,&dummy);
0105         
0106         x=x/equilparam.ap[0];//Mandatory normalization
0107         y=y/equilparam.ap[0];//Mandatory normalization
0108         dxdrho=dxdrho/equilparam.ap[0];//Mandatory normalization
0109         dydrho=dydrho/equilparam.ap[0];//Mandatory normalization
0110         
0111         down[1] = dxdrho*(x-Y[0]) + dydrho*(y - Y[1]);
0112         
0113         ppvald_2D(rayequil.x_fit,rayparam.dS[0],res[0],theta[0],&x,&dxdrho,&dummy,&dummy,&dummy,&dummy);
0114         ppvald_2D(rayequil.y_fit,rayparam.dS[0],res[0],theta[0],&y,&dydrho,&dummy,&dummy,&dummy,&dummy);
0115         
0116         x=x/equilparam.ap[0];//Mandatory normalization
0117         y=y/equilparam.ap[0];//Mandatory normalization
0118         dxdrho=dxdrho/equilparam.ap[0];//Mandatory normalization
0119         dydrho=dydrho/equilparam.ap[0];//Mandatory normalization
0120         
0121         res[1] = dxdrho*(x-Y[0]) + dydrho*(y - Y[1]);
0122         
0123         
0124         save = res[0];
0125         res[0] = res[0] - res[1]*((res[0] - down[0])/(res[1] - down[1]));
0126         down[0] = save;
0127         down[1] = res[1];
0128         
0129         k = k +1;
0130     }
0131     rho[0]=res[0];
0132 }
0133 
0134 void rho1D_newton_analytic(rayequil_format rayequil,rayparam_format rayparam, equilparam_format equilparam, double *rho,double *theta,double *Y)
0135 {
0136     double res[3];
0137 
0138     res[0]=0.05;/*res[0]=rho, res[1] = d(distance^2/2)/drho, res[2]= d2(distance^2/2)drho2*/
0139     double dummy,x,y,dxdrho,dydrho,d2xdrho2,d2ydrho2;
0140     int k = 0;
0141     res[1]=1;
0142     int kmax = 10;
0143     double prec = 1e-9;
0144     while(k < kmax && fabs(res[1]) > prec)
0145     {
0146         ppvald_2D(rayequil.x_fit,rayparam.dS[0],res[0],theta[0],&x,&dxdrho,&dummy,&dummy,&dummy,&d2xdrho2);
0147         ppvald_2D(rayequil.y_fit,rayparam.dS[0],res[0],theta[0],&y,&dydrho,&dummy,&dummy,&dummy,&d2ydrho2);
0148         
0149         x=x/equilparam.ap[0];//Mandatory normalization
0150         y=y/equilparam.ap[0];//Mandatory normalization
0151         dxdrho=dxdrho/equilparam.ap[0];//Mandatory normalization
0152         dydrho=dydrho/equilparam.ap[0];//Mandatory normalization
0153         d2xdrho2=d2xdrho2/equilparam.ap[0];//Mandatory normalization
0154         d2ydrho2=d2ydrho2/equilparam.ap[0];//Mandatory normalization
0155         
0156         res[1] = dxdrho*(x-Y[0]) + dydrho*(y - Y[1]);
0157         
0158         res[2] = d2xdrho2*(x-Y[0]) + d2ydrho2*(y - Y[1]) + pow(dxdrho,2) + pow(dydrho,2);
0159         
0160         res[0] = res[0] - res[1]/res[2];
0161         
0162         k = k + 1;
0163     }
0164     rho[0]=res[0];
0165     mexPrintf("%d\n",k);
0166 }
0167 
0168 void substitution_dichotomy_2D(rayequil_format rayequil,rayparam_format rayparam, equilparam_format equilparam, double *rho,double *theta,double *Y)
0169 {
0170     theta[0] = atan2(Y[1],Y[0]);
0171     double down_t[2],up_t[2],res_t[2];
0172     res_t[0] = theta[0];
0173     down_t[0] = theta[0]*0.9;
0174     up_t[0] = theta[0]*1.1;
0175     double dummy,x,y,dxdtheta,dydtheta,rho_d,rho_m;
0176     int k=0;
0177     int kmax = 60;
0178     double prec = 1e-11;
0179     res_t[1] = 1;
0180     while(k < kmax && fabs(res_t[1]) > prec)
0181     {
0182         k = k + 1;
0183         res_t[0] = (down_t[0] + up_t[0])/2;
0184         rho1D_dichotomy(rayequil,rayparam,equilparam,&rho_d,&down_t[0],Y);
0185         rho1D_dichotomy(rayequil,rayparam,equilparam,&rho_m,&res_t[0],Y);
0186         
0187         ppvald_2D(rayequil.x_fit,rayparam.dS[0],rho_d,down_t[0],&x,&dummy,&dxdtheta,&dummy,&dummy,&dummy);
0188         ppvald_2D(rayequil.y_fit,rayparam.dS[0],rho_d,down_t[0],&y,&dummy,&dydtheta,&dummy,&dummy,&dummy);
0189         x=x/equilparam.ap[0];//Mandatory normalization
0190         y=y/equilparam.ap[0];//Mandatory normalization
0191         dxdtheta=dxdtheta/equilparam.ap[0];//Mandatory normalization
0192         dydtheta=dydtheta/equilparam.ap[0];//Mandatory normalization
0193         down_t[1] = dxdtheta*(x-Y[0]) + dydtheta*(y-Y[1]);
0194         
0195         ppvald_2D(rayequil.x_fit,rayparam.dS[0],rho_m,res_t[0],&x,&dummy,&dxdtheta,&dummy,&dummy,&dummy);
0196         ppvald_2D(rayequil.y_fit,rayparam.dS[0],rho_m,res_t[0],&y,&dummy,&dydtheta,&dummy,&dummy,&dummy);
0197         x=x/equilparam.ap[0];//Mandatory normalization
0198         y=y/equilparam.ap[0];//Mandatory normalization
0199         dxdtheta=dxdtheta/equilparam.ap[0];//Mandatory normalization
0200         dydtheta=dydtheta/equilparam.ap[0];//Mandatory normalization
0201         res_t[1] = dxdtheta*(x-Y[0]) + dydtheta*(y-Y[1]);
0202         
0203         if(down_t[1]*res_t[1]>0)
0204         {
0205             down_t[0] = res_t[0];
0206         }
0207         else
0208         {
0209             up_t[0] = res_t[0];
0210         }
0211     }
0212     theta[0] = res_t[0];
0213     rho1D_dichotomy(rayequil,rayparam,equilparam,rho,theta,Y);
0214 }
0215 
0216 void substitution_newton_numeric_2D(rayequil_format rayequil,rayparam_format rayparam, equilparam_format equilparam, double *rho,double *theta,double *Y)
0217 {
0218     theta[0] = atan2(Y[1],Y[0]);
0219     double down_t[2],res_t[2];
0220     down_t[0] = theta[0];
0221     res_t[0] = 1.00000001*down_t[0];
0222     double dummy,x,y,dxdtheta,dydtheta,rho_d,rho_m;
0223     int kt=0;
0224     int ktmax = 15;
0225     double prec = 1e-9;
0226     res_t[1] = 1;
0227     double save;
0228     
0229     rho1D_newton_analytic(rayequil,rayparam,equilparam,&rho_d,&down_t[0],Y);
0230     ppvald_2D(rayequil.x_fit,rayparam.dS[0],rho_d,down_t[0],&x,&dummy,&dxdtheta,&dummy,&dummy,&dummy);
0231     ppvald_2D(rayequil.y_fit,rayparam.dS[0],rho_d,down_t[0],&y,&dummy,&dydtheta,&dummy,&dummy,&dummy);
0232     x=x/equilparam.ap[0];//Mandatory normalization
0233     y=y/equilparam.ap[0];//Mandatory normalization
0234     dxdtheta=dxdtheta/equilparam.ap[0];//Mandatory normalization
0235     dydtheta=dydtheta/equilparam.ap[0];//Mandatory normalization
0236     down_t[1] = dxdtheta*(x-Y[0]) + dydtheta*(y-Y[1]);
0237     
0238     while(kt < ktmax && fabs(res_t[1]) > prec)
0239     {
0240         rho1D_newton_analytic(rayequil,rayparam,equilparam,&rho_m,&res_t[0],Y);
0241         
0242         ppvald_2D(rayequil.x_fit,rayparam.dS[0],rho_m,res_t[0],&x,&dummy,&dxdtheta,&dummy,&dummy,&dummy);
0243         ppvald_2D(rayequil.y_fit,rayparam.dS[0],rho_m,res_t[0],&y,&dummy,&dydtheta,&dummy,&dummy,&dummy);
0244         x=x/equilparam.ap[0];//Mandatory normalization
0245         y=y/equilparam.ap[0];//Mandatory normalization
0246         dxdtheta=dxdtheta/equilparam.ap[0];//Mandatory normalization
0247         dydtheta=dydtheta/equilparam.ap[0];//Mandatory normalization
0248         
0249         res_t[1] = dxdtheta*(x-Y[0]) + dydtheta*(y-Y[1]);
0250         
0251         save = res_t[0];
0252         res_t[0] = res_t[0] - res_t[1]*((res_t[0] - down_t[0])/(res_t[1] - down_t[1]));
0253         down_t[0] = save;
0254         down_t[1] = res_t[1];
0255         
0256         DUMP(res_t[0]);
0257         
0258         kt = kt + 1;
0259         
0260     }
0261     theta[0] = res_t[0];
0262     rho1D_newton_analytic(rayequil,rayparam,equilparam,rho,theta,Y);
0263 }
0264 
0265 
0266 void substitution_newton_analytic_2D(rayequil_format rayequil,rayparam_format rayparam, equilparam_format equilparam, double *rho,double *theta,double *Y)
0267 {
0268     theta[0] = atan2(Y[1],Y[0]);
0269     double res_t[3];//res_t[0]=theta, res_t[1]=d(distance^2/2)dtheta, res_t[2]=d2(distance^2/2)dtheta2
0270     res_t[0] = theta[0];
0271     double dummy,x,y,dxdtheta,dydtheta,d2xdtheta2,d2ydtheta2,rho_s;
0272     int kt=0;
0273     int ktmax = 15;
0274     double prec = 1e-9;
0275     res_t[1] = 1;
0276     while(kt < ktmax && fabs(res_t[1]) > prec)
0277     {
0278         rho1D_newton_analytic(rayequil,rayparam,equilparam,&rho_s,&res_t[0],Y);
0279         
0280         ppvald_2D(rayequil.x_fit,rayparam.dS[0],rho_s,res_t[0],&x,&dummy,&dxdtheta,&d2xdtheta2,&dummy,&dummy);
0281         ppvald_2D(rayequil.y_fit,rayparam.dS[0],rho_s,res_t[0],&y,&dummy,&dydtheta,&d2ydtheta2,&dummy,&dummy);
0282         
0283         x=x/equilparam.ap[0];//Mandatory normalization
0284         y=y/equilparam.ap[0];//Mandatory normalization
0285         dxdtheta=dxdtheta/equilparam.ap[0];//Mandatory normalization
0286         dydtheta=dydtheta/equilparam.ap[0];//Mandatory normalization
0287         d2xdtheta2=d2xdtheta2/equilparam.ap[0];//Mandatory normalization
0288         d2ydtheta2=d2ydtheta2/equilparam.ap[0];//Mandatory normalization
0289         
0290         res_t[1] = dxdtheta*(x-Y[0]) + dydtheta*(y-Y[1]);
0291         
0292         res_t[2] = d2xdtheta2*(x-Y[0]) + d2ydtheta2*(y-Y[1]) + pow(dxdtheta,2) + pow(dydtheta,2);
0293         
0294         res_t[0] = res_t[0] - res_t[1]/res_t[2];
0295         
0296         DUMP(res_t[0]);
0297         
0298         kt = kt + 1;
0299         
0300     }
0301     theta[0] = res_t[0];
0302     rho1D_newton_analytic(rayequil,rayparam,equilparam,rho,theta,Y);
0303 }
0304 
0305 void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
0306 {
0307     double *x,*y,*rho,*theta,*Y;
0308     
0309     x = mxCalloc(1,sizeof(double));
0310     y = mxCalloc(1,sizeof(double));
0311     Y = mxCalloc(4,sizeof(double));
0312     
0313     /*Declaration of the input parameters*/
0314     
0315     x       =   mxGetPr(x_IN);
0316     y       =   mxGetPr(y_IN);
0317     
0318     /*Declaration of the output parameters*/
0319     
0320     rho_OUT = mxCreateDoubleScalar(mxREAL);
0321     theta_OUT = mxCreateDoubleScalar(mxREAL);
0322     rho = mxGetPr(rho_OUT);
0323     theta = mxGetPr(theta_OUT);
0324     
0325     
0326     // Declaration of the structures
0327     
0328     rayparam_format rayparam;
0329     rayequil_format rayequil;
0330     equilparam_format equilparam;
0331     
0332     // Loading of the magnetic equilibrium
0333     
0334     rayequil.x_fit       = load_grid_2D(equil_IN,"x_fit",0);
0335     rayequil.y_fit       = load_grid_2D(equil_IN,"y_fit",0);
0336     rayparam             = load_ray_param(rayparam_IN);
0337     equilparam           = load_equil_param(equilparam_IN);
0338     
0339     //Scheme
0340     Y[0]=x[0];
0341     Y[1]=y[0];
0342     
0343     substitution_newton_analytic_2D(rayequil,rayparam,equilparam,rho,theta,Y);
0344 
0345     return;
0346 }
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.