fluctstruct_yp

PURPOSE ^

SYNOPSIS ^

function [fluct] = fluctstruct_yp(equil,id,fluct_in)

DESCRIPTION ^

 Calculates the fluctuations structure for C3PO

 INPUT
    - equil: axisymmetric magnetic equilibrium structure
    - id: fluct id 
    - fluct_in: 
       - naequilp: non axisymmetric perturbation of the toroidal MHD equilibrium (dynamic -> ne: density, B: magnetic field), or (static -> magnetic ripple)
       - npar0_lh: perturbation of the launched npar for the LH wave only (simulate possible diffraction of the LH wave in the SOL by edge density fluctuations)

 OUTPUT:

    - fluct: fluctuations or modulations output structure for C3PO

 By Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [fluct] = fluctstruct_yp(equil,id,fluct_in)
0002 %
0003 % Calculates the fluctuations structure for C3PO
0004 %
0005 % INPUT
0006 %    - equil: axisymmetric magnetic equilibrium structure
0007 %    - id: fluct id
0008 %    - fluct_in:
0009 %       - naequilp: non axisymmetric perturbation of the toroidal MHD equilibrium (dynamic -> ne: density, B: magnetic field), or (static -> magnetic ripple)
0010 %       - npar0_lh: perturbation of the launched npar for the LH wave only (simulate possible diffraction of the LH wave in the SOL by edge density fluctuations)
0011 %
0012 % OUTPUT:
0013 %
0014 %    - fluct: fluctuations or modulations output structure for C3PO
0015 %
0016 % By Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr)
0017 %
0018 if nargin <= 2,
0019     fluct = [];
0020     warning('Not enough input arguments in fluctstruct_yp.m: fluctuation structure is empty.');
0021     return;
0022 else
0023     if ~isfield(fluct_in,'naequilp') && ~isfield(fluct_in,'npar0_lh'),
0024         error('Unrecognized fluctuation types: naequilp or npar0_lh.');
0025     end
0026 end
0027 %
0028 % Axisymmetric magnetic equilibrium
0029 %
0030 [equilDKE] = equilibrium_jd(equil,[],0,'linear',1);
0031 %
0032 [xq] = qfactors_dke_yp(1,equilDKE.Rp,equilDKE.ap,equilDKE.Xx,equilDKE.Xy,equilDKE.XBx,equilDKE.XBy,equilDKE.XBphi,equilDKE.xx0,equilDKE.xB0,equilDKE.xBT0,equilDKE.xBp0);%safety factor
0033 if isnan(xq(1)), xq(1) = xq(2);end
0034 %
0035 dtheta = gradient(equilDKE.mtheta);
0036 [dRdtheta,dRdrho] = gradient(equilDKE.Xx + equilDKE.Rp);%XR
0037 dRdtheta = dRdtheta./(ones(length(equilDKE.xrho),1)*dtheta);
0038 dRdrho = dRdrho./(equilDKE.xrho(:)*ones(1,length(dtheta)));
0039 %
0040 [dZdtheta,dZdrho] = gradient(equilDKE.Xy + equilDKE.Zp);%XZ
0041 dZdtheta = dZdtheta./(ones(length(equilDKE.xrho),1)*dtheta);
0042 dZdrho = dZdrho./(equilDKE.xrho(:)*ones(1,length(dtheta)));
0043 %
0044 XR = equilDKE.Xx + equilDKE.Rp;
0045 Xr = sqrt(equilDKE.Xx.^2 + equilDKE.Xy.^2);
0046 Xsalpha = (equilDKE.Xx.*dRdtheta + equilDKE.Xy.*dZdtheta)./Xr./sqrt(dRdtheta.^2 + dZdtheta.^2);
0047 Xsalpha(1,:) = Xsalpha(2,:);
0048 Xsalpha(end,:) = Xsalpha(end-1,:);
0049 Xcalpha = sqrt(1.0 - Xsalpha.^2);
0050 Xcalpha2 = Xcalpha.*Xcalpha;
0051 %
0052 sigmaIp =  sign(equil.psi_apRp(end));% Sign of the plasma current (only valid for tokamak, not for RFP !)
0053 XBP = sqrt(equilDKE.XBx.^2 + equilDKE.XBy.^2);
0054 XB = sqrt(equilDKE.XBx.^2 + equilDKE.XBy.^2 + equilDKE.XBphi.^2);
0055 XP = sigmaIp*XBP./XB;
0056 XP2 = XP.*XP;
0057 XT = equilDKE.XBphi./XB;
0058 XT2 = XT.*XT;
0059 %
0060 %Xrho = equilDKE.xrho(:)*ones(size(equilDKE.mtheta));
0061 %Xtheta = ones(size(equilDKE.xrho(:)))*equilDKE.mtheta;
0062 %
0063 L_rho = equilDKE.ap;
0064 xL_theta = sum(Xr./Xcalpha.*(ones(length(equilDKE.xrho),1)*dtheta),2);
0065 xL_perp = sum(Xr./Xcalpha./abs(XT).*(ones(length(equilDKE.xrho),1)*dtheta),2)./xq(:);
0066 XL_phi = 2*pi*(equilDKE.Xx + equilDKE.Rp);
0067 %
0068 % Build fluct structure for C3PO. One substructure for each fluctuation model
0069 %
0070 fluct.id = id;
0071 fluct.equil_id = equil.id;
0072 %
0073 fluct.rho = equilDKE.xrho(:);
0074 fluct.theta = equilDKE.mtheta;
0075 fluct.phi = [];
0076 %
0077 fluct.xq = xq(:);
0078 %
0079 fluct.L_rho = L_rho;
0080 fluct.xL_theta = xL_theta;
0081 fluct.xL_perp = xL_perp;
0082 fluct.XL_phi = XL_phi;
0083 %
0084 fluct.Xcm = Xr./sqrt(1.0 + (XP2./XT2).*Xcalpha2)/equilDKE.ap;%kpar = 0 fluctuation spectrum model
0085 fluct.Xcm(:,end) = fluct.Xcm(:,1);%for spline-Fourier series expansion in the interpolation procedure
0086 fluct.Xcn = -XR.*(XP./XT).*Xcalpha./sqrt(1.0 + (XP2./XT2).*Xcalpha2)/equilDKE.ap;%kpar = 0 fluctuation spectrum model
0087 fluct.Xcn(:,end) = fluct.Xcn(:,1);%for spline-Fourier series expansion in the interpolation procedure
0088 %
0089 if isfield(fluct_in,'naequilp'), 
0090     %
0091     if isfield(fluct_in.naequilp,'itnmax'),
0092         fluct.itnmax = fluct_in.naequilp.itnmax;% number of iterations
0093     end
0094     %
0095     inefluct_model = 1;
0096     iBfluct_model = 1;
0097     %
0098     if isfield(fluct_in.naequilp,'fluct'),
0099         for i = 1:length(fluct_in.naequilp.fluct.type),        
0100             if strcmp(fluct_in.naequilp.fluct.type{i},'ne'),
0101                 fluct.ne(inefluct_model).model = fluct_in.naequilp.fluct.model(i);%fluctuation model
0102                 %
0103                 fluct.ne(inefluct_model).sigmar_max = fluct_in.naequilp.fluct.sigmar_max(i);%Maximum value of the relative fluctuations variance at the poloidal angle theta = 0 [1,nfluct_types]
0104                 fluct.ne(inefluct_model).sigmar_hwhm = fluct_in.naequilp.fluct.sigmar_hwhm(i);%Radial half width at half maximum of the relative fluctuations variance at the poloidal angle theta = 0 [1,nfluct_types]
0105                 fluct.ne(inefluct_model).sigmar_rho = fluct_in.naequilp.fluct.sigmar_rho(i);%Radial location where the relative fluctuations variance peaks at the poloidal angle theta = 0 [1,nfluct_types]
0106                 fluct.ne(inefluct_model).polmode = fluct_in.naequilp.fluct.polmode(i);%HFS/LFS relative fluctuations variance ratio. No poloidal dependency corresponds to 1 [1,nfluct_types]
0107                 %
0108                 fluct.ne(inefluct_model).epsi_rho = fluct_in.naequilp.fluct.epsi_rho(i);%
0109                 fluct.ne(inefluct_model).epsi_theta = fluct_in.naequilp.fluct.epsi_theta(i);%
0110                 fluct.ne(inefluct_model).epsi_phi = fluct_in.naequilp.fluct.epsi_phi(i);%
0111                 %
0112                 fluct.ne(inefluct_model).lmin = fluct_in.naequilp.fluct.lmin(i);%
0113                 fluct.ne(inefluct_model).mmin = fluct_in.naequilp.fluct.mmin(i);%
0114                 fluct.ne(inefluct_model).nmin = fluct_in.naequilp.fluct.nmin(i);%
0115                 %
0116                 fluct.ne(inefluct_model).lmax = fluct_in.naequilp.fluct.lmax(i);%
0117                 fluct.ne(inefluct_model).mmax = fluct_in.naequilp.fluct.mmax(i);%
0118                 fluct.ne(inefluct_model).nmax = fluct_in.naequilp.fluct.nmax(i);%
0119                 %
0120                 inefluct_model = inefluct_model + 1;
0121                 %
0122             elseif strcmp(fluct_in.naequilp.fluct.type{i},'B'),
0123                 %
0124                 fluct.B(iBfluct_model).model = fluct_in.naequilp.fluct.model(i);%fluctuation model
0125                 %
0126                 fluct.B(iBfluct_model).sigmar_max = fluct_in.naequilp.fluct.sigmar_max(i);%Maximum value of the relative fluctuations variance at the poloidal angle theta = 0 [1,nfluct_types]
0127                 fluct.B(iBfluct_model).sigmar_hwhm = fluct_in.naequilp.fluct.sigmar_hwhm(i);%Radial half width at half maximum of the relative fluctuations variance at the poloidal angle theta = 0 [1,nfluct_types]
0128                 fluct.B(iBfluct_model).sigmar_rho = fluct_in.naequilp.fluct.sigmar_rho(i);%Radial location where the relative fluctuations variance peaks at the poloidal angle theta = 0 [1,nfluct_types]
0129                 fluct.B(iBfluct_model).polmode = fluct_in.naequilp.fluct.polmode(i);%HFS/LFS relative fluctuations variance ratio. No poloidal dependency corresponds to 1 [1,nfluct_types]
0130                 %
0131                 fluct.B(iBfluct_model).epsi_rho = fluct_in.naequilp.fluct.epsi_rho(i);%
0132                 fluct.B(iBfluct_model).epsi_theta = fluct_in.naequilp.fluct.epsi_theta(i);%
0133                 fluct.B(iBfluct_model).epsi_phi = fluct_in.naequilp.fluct.epsi_phi(i);%
0134                 %
0135                 fluct.B(iBfluct_model).lmin = fluct_in.naequilp.fluct.lmin(i);%
0136                 fluct.B(iBfluct_model).mmin= fluct_in.naequilp.fluct.mmin(i);%
0137                 fluct.B(iBfluct_model).nmin = fluct_in.naequilp.fluct.nmin(i);%
0138                 %
0139                 fluct.B(iBfluct_model).lmax = fluct_in.naequilp.fluct.lmax(i);%
0140                 fluct.B(iBfluct_model).mmax= fluct_in.naequilp.fluct.mmax(i);%
0141                 fluct.B(iBfluct_model).nmax = fluct_in.naequilp.fluct.nmax(i);%
0142                 %
0143                 iBfluct_model = iBfluct_model + 1;
0144             end
0145         end
0146     end
0147     %
0148     % Magnetic ripple model
0149     %
0150     if isfield(fluct_in.naequilp,'ripple'), 
0151         %
0152         % B field terms
0153         %
0154         [aBx,aBy,aBz] = bripfield_yp(fluct_in.naequilp.ripple,equilDKE,3);%mode 3 for circular shaped windings cross-section
0155         B = sqrt(equilDKE.XBx.^2 + equilDKE.XBy.^2 + equilDKE.XBphi.^2);%
0156         %
0157         fluct.B(iBfluct_model).model = 0;%Magnetic ripple model
0158         %
0159         fluct.B(iBfluct_model).sigmar_max = [];
0160         fluct.B(iBfluct_model).sigmar_hwhm = [];
0161         fluct.B(iBfluct_model).sigmar_rho = [];
0162         fluct.B(iBfluct_model).polmode = [];
0163         %
0164         fluct.B(iBfluct_model).epsi_rho = [];
0165         fluct.B(iBfluct_model).epsi_theta = [];
0166         fluct.B(iBfluct_model).epsi_phi = [];
0167         %
0168         fluct.B(iBfluct_model).lmin = [];
0169         fluct.B(iBfluct_model).mmin= [];
0170         fluct.B(iBfluct_model).nmin = fluct_in.naequilp.ripple.N;%Number of toroidal magnetic field coils (one harmonic)
0171         %
0172         fluct.B(iBfluct_model).lmax = [];
0173         fluct.B(iBfluct_model).mmax= [];
0174         fluct.B(iBfluct_model).nmax = fluct_in.naequilp.ripple.N;%Number of toroidal magnetic field coils (one harmonic)
0175         %
0176         fluct.B(iBfluct_model).Xfx = aBx./B;
0177         fluct.B(iBfluct_model).Xfy = aBy./B;
0178         fluct.B(iBfluct_model).Xfz = aBz./B;%here z stand for phi the toroidal direction
0179         %
0180         iBfluct_model = iBfluct_model + 1;
0181     end
0182 end
0183 %
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.