CPO2LUKEstruct_yp

PURPOSE ^

LUKE - Function that build a LUKE compliant structure from the corresponding CPO in the ITM MDS+ database

SYNOPSIS ^

function [LUKE_struct] = CPO2LUKEstruct_yp(ver_mds,struct_name,CPO_struct,display_mode)

DESCRIPTION ^

LUKE - Function that build a LUKE compliant structure from the corresponding CPO in the ITM MDS+ database

Function that build a LUKE compliant structure from the corresponding CPO in the ITM MDS+ database

 Input:

   - ver_mds: version of the ITM MDS+ data format (4.09a, 4.10a,...). Important for the conversion to LUKE data format
    - struct_name: CPO name that characterizes its type in the ITM database ('equilibrium', 'coreprof', etc)
    - CPO_struct: CPO data
    - display_mode: display results if any (0: do nothing, > 0: display results)

 Output:

    - LUKE_struct: structure totally or partially compatible with LUKE package (for consistency problems with other ITM CPO's)

by Y.Peysson CEA-IRFM <yves.peysson@cea.fr> and Joan Decker CEA-IRFM (joan.decker@cea.fr)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [LUKE_struct] = CPO2LUKEstruct_yp(ver_mds,struct_name,CPO_struct,display_mode)
0002 %LUKE - Function that build a LUKE compliant structure from the corresponding CPO in the ITM MDS+ database
0003 %
0004 %Function that build a LUKE compliant structure from the corresponding CPO in the ITM MDS+ database
0005 %
0006 % Input:
0007 %
0008 %   - ver_mds: version of the ITM MDS+ data format (4.09a, 4.10a,...). Important for the conversion to LUKE data format
0009 %    - struct_name: CPO name that characterizes its type in the ITM database ('equilibrium', 'coreprof', etc)
0010 %    - CPO_struct: CPO data
0011 %    - display_mode: display results if any (0: do nothing, > 0: display results)
0012 %
0013 % Output:
0014 %
0015 %    - LUKE_struct: structure totally or partially compatible with LUKE package (for consistency problems with other ITM CPO's)
0016 %
0017 %by Y.Peysson CEA-IRFM <yves.peysson@cea.fr> and Joan Decker CEA-IRFM (joan.decker@cea.fr)
0018 %
0019 [dummy,dummy,dummy,dummy,dummy,dummy,dummy,dummy,clum] = pc_dke_yp;%Speed of light (m/s)
0020 %
0021 if nargin < 3,error('Not enough input arguments in CPO2LUKE_struct_yp.m');end
0022 %
0023 if nargin == 3,
0024     display_mode = 0;%display option
0025 end
0026 %
0027 p_opt = -1;%printing option
0028 %
0029 mds_dataversion0 = str2num(ver_mds(1:end-1));%for using numbers
0030 mds_dataversion1 = ver_mds(end:end);%for using
0031 %
0032 if strcmp('equilibrium',struct_name),
0033     %
0034     if mds_dataversion0 >= 4.10,
0035         cocos = CPO_struct.datainfo.cocos;
0036     else
0037         cocos = '';
0038     end
0039     %
0040     grid_type = CPO_struct.profiles_2d.grid_type;%CPO is a java structure and it must be transfered to Matlab data first before any operation
0041     %
0042     if round(str2double(grid_type(1))) == 1,%rectangular grid
0043         xyBR = CPO_struct.profiles_2d.br;
0044         xyBZ = CPO_struct.profiles_2d.bz;
0045         xyBPHI = CPO_struct.profiles_2d.bphi;    
0046         %
0047         xypsi = CPO_struct.profiles_2d.psi;
0048         xypsi = xypsi/2/pi;
0049         psi = CPO_struct.profiles_1d.psi;
0050         psi = psi/2/pi;%Poloidal flux (with 1/2/pi)
0051         psia = CPO_struct.global_param.psi_bound;
0052         psia = psia/2/pi;
0053         psi0 = CPO_struct.global_param.psi_ax;
0054         psi0 = psi0/2/pi;
0055         %
0056         Rp = CPO_struct.global_param.mag_axis.position.r;    
0057         Zp = CPO_struct.global_param.mag_axis.position.z;
0058         %
0059         x = CPO_struct.profiles_2d.grid.dim1;
0060         x = x - Rp;
0061         y = CPO_struct.profiles_2d.grid.dim2;
0062         y = y - Zp;
0063         %
0064         Btor0 = CPO_struct.global_param.toroid_field.b0;%vacuum toroidal field (T) (for the sign only)
0065         %
0066         Ip = CPO_struct.global_param.i_plasma;%plasma current (A) (for the sign only)
0067         qspi = CPO_struct.profiles_1d.q;%signed safety factor
0068         %
0069         if psi0 < psia, %IPHI is positive, PHI clockwise from top
0070             sigma = 1;
0071         else %IPHI is negative, PHI clockwise from top
0072             sigma = -1;
0073         end
0074         %
0075         xypsi = xypsi - psi0;%psi = 0 on the magnetic axis (arbitrary choice)
0076         psi = psi - psi0;
0077         psia = psia - psi0;
0078         %
0079         xypsin = xypsi/psia;%normalized poloidal magnetic flux
0080         psin = psi(:)'/psia;
0081         npsin = length(psin);%number of radial points (normalized grid point)
0082         %
0083         ntheta = 65;%number of theta points (poloidal angle)
0084         opt.sfac = 24;%have to increased when ray goes close to the magnetic axis
0085         opt.axenforce = true;
0086         %
0087         % (R,Z) to (psi,theta) grid conversion (two methods)
0088         %
0089         if isempty(xyBR) || isempty(xyBZ),%calculations of the poloidal B field from the poloidal magnetic flux
0090             %
0091             equil_magnetic_data.x = x;
0092             equil_magnetic_data.y = y;
0093             equil_magnetic_data.Rp = Rp;
0094             equil_magnetic_data.Zp = Zp;
0095             equil_magnetic_data.xypsin = xypsin;
0096             equil_magnetic_data.psia = psia;
0097             equil_magnetic_data.xBt = -xyBPHI(:,1);%sign change -> ITM: BPHI > 0 if direction counter clockwise, LUKE: BPHI > 0 direction is clockwise (both view from the top)
0098             equil_magnetic_data.psin = psin;
0099             equil_magnetic_data.pq = -qspi;%sign change -> ITM: Ip > 0 if direction counter clockwise, LUKE: Ip > 0 direction is clockwise (both view from the top)
0100             %
0101             [LUKE_struct,prho,pspsin,pspsinT] = equil_magnetic_rz2pt_jd(equil_magnetic_data,display_mode,p_opt,'','',npsin,ntheta,opt);
0102             %
0103         else,%direct interpolation of the magnetic field
0104             interp_method = 'spline';%interpolation method for grid conversion.
0105             %
0106             [ptx,pty,ptdpsindx,ptdpsindy,psin,theta] = rz2pt_jd(x,y,xypsin,psin,ntheta,interp_method,display_mode,p_opt,'',opt);
0107             %
0108             ptBPHI = -griddata(x,y,xyBPHI',ptx,pty);%sign change -> ITM: BPHI > 0 if direction counter clockwise, LUKE: BPHI > 0 direction is clockwise (both view from the top)
0109             %
0110             ap = ptx(npsin,1);%Plasma minor radius on the LFS midplane
0111             psi_apRp = psin*psia*ap/Rp;%The poloidal flux coordinate is normalized by the aspect ratio
0112             ptR = Rp + ptx;
0113             %
0114             % Note that in LUKE we have B = F.gradphi + gradphi x gradpsi.
0115             %
0116             ptBx = griddata(x,y,xyBR',ptx,pty);
0117             ptBy = griddata(x,y,xyBZ',ptx,pty);
0118             %
0119             ptBP = sqrt(ptBx.^2 +ptBy.^2);
0120             ptB = sqrt(ptBP.^2 + ptBPHI.^2);
0121             %
0122             prho = ptx(:,1).'/ap;
0123             %
0124             if display_mode == 2, 
0125                 %
0126                 rx = [-fliplr(prho(2:npsin)),prho];
0127                 rBx = [fliplr(ptBx(2:npsin,(ntheta+1)/2)'),ptBx(:,1)'];
0128                 rBy = [fliplr(ptBy(2:npsin,(ntheta+1)/2)'),ptBy(:,1)'];
0129                 rBPHI = [fliplr(ptBPHI(2:npsin,(ntheta+1)/2)'),ptBPHI(:,1)'];
0130                 rBP = sqrt(rBx.^2 + rBy.^2);
0131                 rB = sqrt(rBP.^2 + rBPHI.^2);
0132                 %
0133                 style = '-';
0134                 marker = 'none';
0135                 color1 = 'r';
0136                 color2 = 'b';
0137                 color3 = 'k';
0138                 width = 2;
0139                 siz = 20;
0140                 red = 0.9;
0141                 lspace = 0.7;
0142                 bspace = 0.6;
0143                 %
0144                 xlab = 'r/a';
0145                 ylab = 'B(T)';
0146                 xlim = [-1,1];
0147                 ylim = NaN;
0148                 leg = ['B_T';'B_P';'B  '];
0149                 %
0150                 figure_yp('','B field radial profiles','clf');
0151                 graph1D_jd(rx,abs(rBPHI),0,0,'','','',NaN,xlim,ylim,style,marker,color1,width,siz,gca,red,lspace,bspace);
0152                 graph1D_jd(rx,rBP,0,0,'','','',NaN,xlim,ylim,style,marker,color2,width);
0153                 graph1D_jd(rx,rB,0,0,xlab,ylab,'',leg,xlim,ylim,style,marker,color3,width);
0154                 %
0155                 figure_yp('','B field poloidal dependence','clf');
0156                 graph1D_jd(theta,ptB(2,:),0,0,'','','',NaN,'','',style,marker,'b',width);hold on
0157                 plot(theta,ptB([2:end],:))
0158                 legend(['r/a = ',num2str(prho(2))])
0159                 axis([0,2*pi,0,ceil(max(max(ptB([2:end],:)))*10)/10])
0160                 xlabel('\theta (rd)'),ylabel('B(T)')
0161             end
0162             %
0163             ptx(:,end) = ptx(:,1);
0164             pty(:,end) = pty(:,1);
0165             ptBx(:,end) = ptBx(:,1);
0166             ptBy(:,end) = ptBy(:,1);
0167             ptBPHI(:,end) = ptBPHI(:,1);
0168             %
0169             % equil structure
0170             %
0171             LUKE_struct.psi_apRp = psi_apRp;
0172             LUKE_struct.theta = theta;
0173             LUKE_struct.ptx = ptx;
0174             LUKE_struct.pty = pty;
0175             LUKE_struct.ptBx = ptBx;
0176             LUKE_struct.ptBy = ptBy;
0177             LUKE_struct.ptBPHI = ptBPHI;
0178             LUKE_struct.Rp = Rp;
0179             LUKE_struct.Zp = Zp;
0180         end
0181         %
0182         % For writing results in CPO in a consistent manner with ITM definitions
0183         %
0184         LUKE_struct.ITM.psi0 = psi0;
0185         LUKE_struct.ITM.psia = psia;
0186         LUKE_struct.ITM.sign_BHI = sign(Btor0);
0187         LUKE_struct.ITM.sign_Ip = sign(Ip);
0188     else
0189         error('This grid is not yet recognized in CPO2LUKE_struct_yp.m.')
0190     end
0191     %
0192     LUKE_struct.id = CPO_struct.datainfo.whatref.machine;
0193     %
0194 elseif strcmp('coreprof',struct_name),
0195     %
0196     LUKE_struct.psi = CPO_struct.psi.value;% Unfortunately, the poloidal psi for the 1-D profiles is not consistent with the psi defined in the equil structure (number of grid points, min and max !!!)... who named that CPO's, Consistent Physical Object....
0197     LUKE_struct.psi = LUKE_struct.psi(:)';% Unfortunately, the poloidal psi for the 1-D profiles is not consistent with the psi defined in the equil structure (number of grid points, min and max !!!)... who named that CPO's, Consistent Physical Object....
0198     LUKE_struct.pTe = CPO_struct.te.value;%keV units
0199     LUKE_struct.pTe = LUKE_struct.pTe(:)'/1000;%keV units
0200     LUKE_struct.pne = CPO_struct.ne.value;
0201     LUKE_struct.pne = LUKE_struct.pne(:)';
0202     LUKE_struct.pzTi = CPO_struct.ti.value;%keV units
0203     LUKE_struct.pzTi = transpose(LUKE_struct.pzTi)/1000;%keV units
0204     LUKE_struct.pzni = CPO_struct.ni.value;
0205     LUKE_struct.pzni = transpose(LUKE_struct.pzni);
0206     LUKE_struct.zZi = CPO_struct.composition.zn;
0207     LUKE_struct.zZi = LUKE_struct.zZi(:)';
0208     LUKE_struct.zmi = CPO_struct.composition.amn;
0209     LUKE_struct.zmi = LUKE_struct.zmi(:)';
0210     %
0211     if isempty(LUKE_struct.zZi) || isempty(LUKE_struct.zmi) || isempty(LUKE_struct.pzTi) || isempty(LUKE_struct.pzni),
0212         %
0213         info_dke_yp(6,'Ion parameters are not consistently described in the ITM MDS+ database. 50-50 deuterium-tritium with electroneutrality is enforced.')
0214         %
0215         LUKE_struct.zZi = [1,1,1]; 
0216         LUKE_struct.zmi = [1,2,3];
0217         %
0218         if ~isempty(LUKE_struct.pTe) && ~isempty(LUKE_struct.pne),
0219             LUKE_struct.pzTi = ones(3,1)*LUKE_struct.pTe;
0220             LUKE_struct.pzni = ones(3,1)*LUKE_struct.pne;
0221             LUKE_struct.pzni(1,:) = 0;
0222             LUKE_struct.pzni(2,:) = LUKE_struct.pzni(2,:)/2;
0223             LUKE_struct.pzni(3,:) = LUKE_struct.pzni(3,:)/2;
0224         else
0225             info_dke_yp(6,'Both ion and electrons parameters are not consistently described. Please check the ITM MDS+ database.')
0226         end
0227     end
0228     %
0229 elseif strcmp('antennas',struct_name),
0230     %
0231     if mds_dataversion0 == 4.09,
0232         %
0233         for ia = 1:length(CPO_struct.antenna_unit),
0234             if ~isempty(CPO_struct.antenna_unit(ia).antenna_ec.name),
0235                 LUKE_struct{ia}.id = CPO_struct.antenna_unit(ia).antenna_ec.name;
0236                 LUKE_struct{ia}.type = 'EC';
0237                 LUKE_struct{ia}.omega_rf = CPO_struct.antenna_unit(ia).antenna_ec.frequency;%wave angular frequency (rad/s)
0238                 LUKE_struct{ia}.omega_rf = LUKE_struct{ia}.omega_rf*2*pi;
0239                 %
0240                 LUKE_struct{ia}.yR_L = CPO_struct.antenna_unit(ia).antenna_ec.position.r;
0241                 LUKE_struct{ia}.yZ_L = CPO_struct.antenna_unit(ia).antenna_ec.position.z;
0242                 LUKE_struct{ia}.yphi_L = CPO_struct.antenna_unit(ia).antenna_ec.position.phi;    
0243                 %
0244                 % Launch angles definitions for ITM are identical to those for Prater benchmark (see ITER simulation example in LUKE)
0245                 %
0246                 beta = CPO_struct.antenna_unit(ia).antenna_ec.launchangles.beta;
0247                 alpha = CPO_struct.antenna_unit(ia).antenna_ec.launchangles.alpha;
0248                 LUKE_struct{ia}.yalpha_L = sign(beta)*pi - atan(tan(beta)/cos(alpha));
0249                 LUKE_struct{ia}.ybeta_L = pi - acos(cos(beta)*sin(alpha));            
0250                 %
0251                 LUKE_struct{ia}.yP_L = CPO_struct.antenna_unit(ia).antenna_ec.power.value;
0252                 LUKE_struct{ia}.dNpar0 = NaN;
0253                 LUKE_struct{ia}.ns = 1000;
0254                 LUKE_struct{ia}.method1 = 'pchip';%1D interpolation
0255                 LUKE_struct{ia}.method2 = 'cubic';%2D interpolation
0256                 %
0257                 frequency = CPO_struct.antenna_unit(ia).antenna_ec.frequency;
0258                 alpha = frequency*pi/clum;%for the Rayleigh length
0259                 waist = CPO_struct.antenna_unit(ia).antenna_ec.beam.spot.waist;
0260                 invcurvrad = CPO_struct.antenna_unit(ia).antenna_ec.beam.phaseellipse.invcurvrad;
0261                 %
0262                 yw = alpha^2*waist(1)^4*invcurvrad(1);%WARNING: circular shape only.
0263                 %
0264                 LUKE_struct{ia}.w0 = waist(1)/sqrt(1 + yw*invcurvrad(1));
0265                 LUKE_struct{ia}.z_L = -yw/(1 + yw*invcurvrad(1));                        
0266                 %
0267                 LUKE_struct{ia}.mmode = -CPO_struct.antenna_unit(ia).antenna_ec.mode;%LUKE definitions -> O/X mode: -1/+1. ITM ones are opposite
0268             end
0269             %
0270             if ~isempty(CPO_struct.antenna_unit(ia).antenna_lh.name),
0271                 disp('WARNING: LH antenna not yet considered in ITM')
0272             end
0273         end
0274     elseif mds_dataversion0 >= 4.10,
0275        %
0276         if isfield(CPO_struct,'antenna_ec'),
0277             for ia = 1:length(CPO_struct.antenna_ec),
0278                 if ~isempty(CPO_struct.antenna_ec(ia).name),
0279                     LUKE_struct{ia}.id = CPO_struct.antenna_ec(ia).name;
0280                     LUKE_struct{ia}.type = 'EC';
0281                     frequency = CPO_struct.antenna_ec(ia).frequency;
0282                     LUKE_struct{ia}.omega_rf = frequency*2*pi;%wave angular frequency (rad/s)
0283                     %
0284                     LUKE_struct{ia}.yR_L = CPO_struct.antenna_ec(ia).position.r;
0285                     LUKE_struct{ia}.yZ_L = CPO_struct.antenna_ec(ia).position.z;
0286                     LUKE_struct{ia}.yphi_L = CPO_struct.antenna_ec(ia).position.phi;    
0287                     %
0288                     % Launch angles definitions for ITM are identical to those for Prater benchmark (see ITER simulation example in LUKE)
0289                     %
0290                     alpha = CPO_struct.antenna_ec(ia).launchangles.alpha;
0291                     beta = CPO_struct.antenna_ec(ia).launchangles.beta;
0292                     %
0293                     LUKE_struct{ia}.yalpha_L = sign(beta)*pi - atan(tan(beta)/cos(alpha));
0294                     LUKE_struct{ia}.ybeta_L = pi - acos(cos(beta)*sin(alpha));            
0295                     %
0296                     LUKE_struct{ia}.yP_L = CPO_struct.antenna_ec(ia).power.value;
0297                     LUKE_struct{ia}.dNpar0 = NaN;
0298                     LUKE_struct{ia}.ns = 1000;
0299                     LUKE_struct{ia}.method1 = 'pchip';%1D interpolation
0300                     LUKE_struct{ia}.method2 = 'cubic';%2D interpolation
0301                     %
0302                     frequency = CPO_struct.antenna_ec(ia).frequency;
0303                     alpha = frequency*pi/clum;%for the Rayleigh length
0304                     %
0305                     size = CPO_struct.antenna_ec(ia).beam.spot.size;
0306                     invcurvrad = CPO_struct.antenna_ec(ia).beam.phaseellipse.invcurvrad;
0307                     %
0308                     yw = alpha^2*size(1)^4*invcurvrad(1);%WARNING: circular shape only.
0309                     %
0310                     LUKE_struct{ia}.w0 = size(1)/sqrt(1 + yw*invcurvrad(1));
0311                     LUKE_struct{ia}.z_L = -yw/(1 + yw*invcurvrad(1));                        
0312                     %
0313                     LUKE_struct{ia}.mmode = -CPO_struct.antenna_ec(ia).mode;%LUKE definitions -> O/X mode: -1/+1. ITM ones are opposite
0314                 else
0315                     
0316                 end
0317             end
0318         else
0319             disp('There is no EC antenna in the MDS+ database.')
0320         end
0321         %
0322          if isfield(CPO_struct,'antenna_lh'),
0323             for ia = 1:length(CPO_struct.antenna_lh),
0324                 if ~isempty(CPO_struct.antenna_lh(ia).name),
0325                     %
0326                 else
0327                     
0328                 end
0329             end
0330         else
0331             disp('There is no LH antenna in the MDS+ database.')
0332          end
0333         %
0334         if isfield(CPO_struct,'antenna_ic'),
0335             for ia = 1:length(CPO_struct.antenna_ic),
0336                 if ~isempty(CPO_struct.antenna_ic(ia).name),
0337                     %
0338                 else
0339                     
0340                 end
0341             end
0342         else
0343             disp('There is no IC antenna in the MDS+ database.')
0344         end    
0345     else
0346         error(['The version ',ver_mds,' of data in EUITM MDS+ database is not recognized in CPO2LUKEstruct_yp.m'])
0347     end
0348 else
0349     error(['There is no CPO with the name ',struct_name,' in the ITM MDS+ database']);
0350 end
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.