qfactors_dke_yp

 Calculation of the different q factors needed for solving the bounce averaged drift-kinetic-equation for
 arbitrary equilibrium (see the notice for explanations concerning notations)

   INPUT:

       - Plasma equilibrium parameters 

   OUTPUT:
   
       - xq (safety factor) [1,nr]
       - xqtilde [1,nr]
       - xqbar [1,nr]
       - xqhat [1,nr]
       - xqcronos [1,nr]

 by Yves Peysson <yves.peysson@cea.fr> (CEA/DRFC) and Joan Decker <jodecker@alum.mit.edu> (MIT/RLE)

0001 function [xq,xqtilde,xqbar,xqhat,xqcronos,xB0] = qfactors_dke_yp(equil_mode,Rp,ap,Xx,Xy,XBx,XBy,XBphi,xx0,xB0,xBT0,xBp0)
0002 %
0003 % Calculation of the different q factors needed for solving the bounce averaged drift-kinetic-equation for
0004 % arbitrary equilibrium (see the notice for explanations concerning notations)
0005 %
0006 %   INPUT:
0007 %
0008 %       - Plasma equilibrium parameters
0009 %
0010 %   OUTPUT:
0011 %
0012 %       - xq (safety factor) [1,nr]
0013 %       - xqtilde [1,nr]
0014 %       - xqbar [1,nr]
0015 %       - xqhat [1,nr]
0016 %       - xqcronos [1,nr]
0017 %
0018 % by Yves Peysson <yves.peysson@cea.fr> (CEA/DRFC) and Joan Decker <jodecker@alum.mit.edu> (MIT/RLE)
0019 %
0020 if equil_mode == 1,
0021     Xx = Xx.';
0022     Xy = Xy.';
0023     XBx = XBx.';
0024     XBy = XBy.';
0025     XBphi = XBphi.';
0026     %
0027     [nt,np] = size(Xx);
0028     %
0029     theta = linspace(0,2*pi,nt); % Angular grid
0030     Xtheta = repmat(theta',[1,np]);
0031     %
0032     XBT = abs(XBphi); % Toroidal B-field
0033     XBP = sqrt(XBx.^2 + XBy.^2); % Poloidal B-field
0034     XB = sqrt(XBT.^2 + XBP.^2); % Total B-field
0035     %
0036     % Half grid in theta for integration
0037     %
0038     if any(isinf(XB(:))),% in case of inf, the half grid is taken as every other point
0039         XMBx = XBx(2:2:nt-1,:);
0040         XMBy = XBy(2:2:nt-1,:);
0041         XMBphi = XBphi(2:2:nt-1,:);
0042         XMB = XB(2:2:nt-1,:);
0043         %
0044         xB0 = min(XMB);
0045         %
0046         XMBmin = ones((nt-1)/2,1)*xB0;
0047         %
0048         XMx = Xx(2:2:nt-1,:);
0049         XMy = Xy(2:2:nt-1,:);
0050         %
0051         XMdt = diff(Xtheta(1:2:end,:));
0052     else
0053         XMBx = (XBx(1:nt-1,:) + XBx(2:nt,:))/2;
0054         XMBy = (XBy(1:nt-1,:) + XBy(2:nt,:))/2;
0055         XMBphi = (XBphi(1:nt-1,:) + XBphi(2:nt,:))/2;
0056         XMB = sqrt(XMBx.^2 + XMBy.^2 + XMBphi.^2);
0057         %
0058         xB0 = min(XMB);
0059         %
0060         XMBmin = ones(nt-1,1)*xB0;
0061         %
0062         XMx = (Xx(1:nt-1,:) + Xx(2:nt,:))/2;
0063         XMy = (Xy(1:nt-1,:) + Xy(2:nt,:))/2;
0064         %
0065         XMdt = diff(Xtheta);
0066     end
0067     %
0068     XMr2 = XMx.^2 + XMy.^2;
0069     XMrBpcos = abs(XMBx.*XMy - XMBy.*XMx);
0070     %
0071     %q (safety factor)
0072     %
0073     XMF = abs(XMBphi).*XMr2./XMrBpcos./(1 + XMx/Rp)/ap;
0074     xq = squeeze(sum(XMF.*XMdt))/(2*pi);
0075     %
0076     %qtilde
0077     %
0078     XMF = XMB.*XMr2./XMrBpcos./ap;
0079     xqtilde = squeeze(sum(XMF.*XMdt))/(2*pi);
0080     %
0081     %qbar
0082     %
0083     XMF = XMBmin.*XMr2./XMrBpcos./(1 + XMx/Rp)/ap;
0084     xqbar = squeeze(sum(XMF.*XMdt))/(2*pi);
0085     %
0086     %qhat
0087     %
0088     XMF = XMBmin.*XMr2./XMrBpcos./ap;
0089     xqhat = squeeze(sum(XMF.*XMdt))/(2*pi);
0090     %
0091     %qcronos
0092     %
0093     XMF = XMB.*XMB.*XMr2./XMrBpcos./XMBmin./ap;
0094     xqcronos = squeeze(sum(XMF.*XMdt))/(2*pi);
0095     %
0096 else
0097     xrho = xx0/ap; %note: for the sake of bounce-averaged calculations in the
0098                           %      analytical limit of circ. conc.
0099                           %      flux-surfaces, we need theta0 = 0, thus
0100                           %      this definition of xrho. Of course, when
0101                           %      the configuration is circ. conc. ro start
0102                           %      with, the two definitions coincide.
0103     xia = xrho*ap/Rp;
0104     xq = xrho.*xBT0./xBp0./sqrt(1 - xia.^2);
0105     xqtilde = xrho.*xB0./xBp0;
0106     xqbar = xrho.*xB0./xBp0./(1 + xia);
0107     xqhat = xqbar;
0108     xqcronos = xrho.*xB0.*sqrt((1+xia)./(1-xia))./xBp0;
0109 end
0110 %
0111