xy2rhotheta_gb

 Find the curvilinear couple associated to a cartesian one

 INPUT
    - equil_fit: interpolated structure
    - rho : Normalized radial coordinate
    - theta : Poloidal angle
    - kmax : Number of iterations of the algorithm

0001 function [rho,theta,rho_0,theta_0] = xy2rhotheta_gb(equil_fit,x,y)
0002 %
0003 % Find the curvilinear couple associated to a cartesian one
0004 %
0005 % INPUT
0006 %    - equil_fit: interpolated structure
0007 %    - rho : Normalized radial coordinate
0008 %    - theta : Poloidal angle
0009 %    - kmax : Number of iterations of the algorithm
0010 
0011 % OUTPUT:
0012 %    - x : Cartesian coordinate x
0013 %    - y : Cartesian coordinate y
0014 
0015 % By Guillaume Brochard (CEA-IRFM, guillaume.brochard@cea.fr)
0016 %
0017 %
0018 function [rho_test] = rho1D(theta_test)%
0019    down = [0,0];
0020    up = [1,0];
0021    res = [0,1];
0022    k=0;
0023    kmax = 60;%Maximum number of iterations allowed
0024    prec = 1e-10;%Precision required
0025    while(k<kmax && abs(res(2))>prec)
0026       k = k+1;
0027       data1 = equilval_yp(equil_fit,down(1),theta_test);
0028       data2 = equilval_yp(equil_fit,up(1),theta_test);
0029       down(2) = data1.dxdrho*(data1.x - x) + data1.dydrho*(data1.y - y);
0030       up(2) = data2.dxdrho*(data2.x - x) + data2.dydrho*(data2.y - y);
0031       res(1) = (down(1) + up(1))/2;
0032       data3 = equilval_yp(equil_fit,res(1),theta_test);
0033       res(2) = data3.dxdrho*(data3.x - x) + data3.dydrho*(data3.y - y);
0034       if(down(2)*res(2)>0)
0035          down = res;
0036       else
0037          up = res;
0038       end
0039    end
0040    rho_test = res(1);
0041 end
0042 %
0043 %
0044 %
0045 theta_0 = atan2(y,x);
0046 rho_0 = rho1D(theta_0);
0047 down_t = [theta_0*(9/10),0];
0048 up_t = [theta_0*(11/10),0];
0049 res_t = [0,1];
0050 l=0;
0051 %
0052 prec = 1e-11;%Precision required
0053 lmax = 60;%Maximum number of iterations allowed
0054 while(l < lmax && abs(res_t(2)) > prec)
0055    l=l+1;
0056    res_t(1) = (down_t(1) + up_t(1))/2;
0057    rho_d = rho1D(down_t(1));
0058    rho_m = rho1D(res_t(1));
0059    data1_t = equilval_yp(equil_fit,rho_d,down_t(1));
0060    data2_t = equilval_yp(equil_fit,rho_m,res_t(1));
0061    down_t(2)= data1_t.dxdtheta*(data1_t.x - x) + data1_t.dydtheta*(data1_t.y - y);
0062    res_t(2) = data2_t.dxdtheta*(data2_t.x - x) + data2_t.dydtheta*(data2_t.y - y);
0063    if(down_t(2)*res_t(2)>0)
0064        down_t = res_t;
0065    else
0066        up_t = res_t;
0067    end
0068 end
0069 %
0070 theta = res_t(1);
0071 rho = rho1D(theta);
0072 %
0073 end