bremchord_dke_yp

 R5-X2 - calculation of the chord characteristics for bremsstrahlung tomography in a torus along unit vector g from LUKE grid structure

 This function calculates chord positions for bremsstrahlung tomography in a torus along unit vector g from LUKE grid structure

 INPUTS

       - equilHXR: magnetic equilibrium structure optimized for HXR calculations
       - radialHXR: Radial grid structure used by LUKE
       - hxr: Bremsstrahlung camera structure
       - hxrparam : Bremsstrahlung calculation and display parameter structure

 OUTPUTS
       
       - Zbremchord: data structure for the bremsstrahlung diagnostic (for display principaly)
          * Zbremchord.vschords: chords parameter [16,n,nc]
          * Zbremchord.psi_hxr: magnetic flux surface coordinate where the distribution function is calculated for the bremsstrahlung emission [npsi,ntheta]
          * Zbremchord.theta_hxr: poloidal angle values where the distribution function is calculated for the bremsstrahlung emission [npsi,ntheta]
          * Zbremchord.PSI_hxr: B/Bmin where the distribution function is calculated for the bremsstrahlung emission [npsi,ntheta]
          * Zbremchord.mask_hxr: locations where the distribution function is calculated for the bremsstrahlung emission [npsi,ntheta]

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

WARNING: the definitions used in the DKE code, and here crucial for HXR calculations are the following
         1) the toroidal frame (R,Z,phi) is direct (phi is positive clockwise when the tokamak is seen from the top)
         2) the poloidal frame (r,theta,phi) is direct
         3) (psi,s,phi) must be direct : if Ip.phi > 0, then B.theta > 0 and psi > 0, and s.theta > 0 (s as the same direction of theta)
         4) the sign of psi is the sign of Ip.phi
 With the definitions above, if Ip is positive, psi is turned outward, Btheta is counterclockwise poloidally, and s too. The Ampere theorem is OK !
 The HXR diagnostic arrangement is taking this frame of reference.

0001 function Zbremchord = bremchord_dke_yp(equilHXR,radialHXR,hxr,hxrparam)
0002 %
0003 % R5-X2 - calculation of the chord characteristics for bremsstrahlung tomography in a torus along unit vector g from LUKE grid structure
0004 %
0005 % This function calculates chord positions for bremsstrahlung tomography in a torus along unit vector g from LUKE grid structure
0006 %
0007 % INPUTS
0008 %
0009 %       - equilHXR: magnetic equilibrium structure optimized for HXR calculations
0010 %       - radialHXR: Radial grid structure used by LUKE
0011 %       - hxr: Bremsstrahlung camera structure
0012 %       - hxrparam : Bremsstrahlung calculation and display parameter structure
0013 %
0014 % OUTPUTS
0015 %
0016 %       - Zbremchord: data structure for the bremsstrahlung diagnostic (for display principaly)
0017 %          * Zbremchord.vschords: chords parameter [16,n,nc]
0018 %          * Zbremchord.psi_hxr: magnetic flux surface coordinate where the distribution function is calculated for the bremsstrahlung emission [npsi,ntheta]
0019 %          * Zbremchord.theta_hxr: poloidal angle values where the distribution function is calculated for the bremsstrahlung emission [npsi,ntheta]
0020 %          * Zbremchord.PSI_hxr: B/Bmin where the distribution function is calculated for the bremsstrahlung emission [npsi,ntheta]
0021 %          * Zbremchord.mask_hxr: locations where the distribution function is calculated for the bremsstrahlung emission [npsi,ntheta]
0022 %
0023 % by Y. Peysson (CEA/IRFM) <yves.peysson@cea.fr> and J. Decker (CEA/IRFM) <joan.decker@cea.fr>
0024 %
0025 %WARNING: the definitions used in the DKE code, and here crucial for HXR calculations are the following
0026 %         1) the toroidal frame (R,Z,phi) is direct (phi is positive clockwise when the tokamak is seen from the top)
0027 %         2) the poloidal frame (r,theta,phi) is direct
0028 %         3) (psi,s,phi) must be direct : if Ip.phi > 0, then B.theta > 0 and psi > 0, and s.theta > 0 (s as the same direction of theta)
0029 %         4) the sign of psi is the sign of Ip.phi
0030 % With the definitions above, if Ip is positive, psi is turned outward, Btheta is counterclockwise poloidally, and s too. The Ampere theorem is OK !
0031 % The HXR diagnostic arrangement is taking this frame of reference.
0032 %
0033 if nargin < 4,
0034    error('Not enough input arguments in bremchord_dke_yp')
0035 end
0036 %
0037 method = hxrparam.method;
0038 ns = hxrparam.ns;
0039 %
0040 xpsin_S = radialHXR.xpsin_S_dke;
0041 %
0042 %   This function interpolates the equilibrium data on the radial grid.
0043 %
0044 Rp = equilHXR.Rp;%Plasma major radius (m)
0045 Zp = equilHXR.Zp;%Plasma vertical position (m)
0046 ap = equilHXR.ap;%Plasma minor radius (m)
0047 theta = equilHXR.mtheta;
0048 ppsin = equilHXR.xpsin;
0049 ptR = equilHXR.Xx + Rp;
0050 ptZ = equilHXR.Xy + Zp;
0051 ptBR = equilHXR.XBx;
0052 ptBZ = equilHXR.XBy;
0053 ptBphi = equilHXR.XBphi;
0054 ptBT = ptBphi;
0055 ptBp = sqrt(ptBR.^2 + ptBZ.^2);
0056 ptB = sqrt(ptBR.^2 + ptBZ.^2 + ptBphi.^2);
0057 ptB0 = ptB(:,1);
0058 ptR0 = equilHXR.Xx(:,1) + Rp;
0059 ptBT0 = ptBT(:,1);
0060 ptBp0 = ptBp(:,1);
0061 ptrho = equilHXR.xrho;
0062 %
0063 nr = length(ppsin);
0064 nr_S = length(xpsin_S);
0065 nt = length(theta);
0066 %
0067 if nr == 1,
0068     error('Bremsstrahlung calculations cannot be performed when LUKE results are obtained for a single radial position !');
0069 end
0070 %
0071 ns = abs(ns);%Force ns to be always positive
0072 nc = length(hxr.alpha_hxr);%chord number
0073 %
0074 xq = ap*qfactors_dke_yp(1,Rp,ap,ptR-Rp,ptZ-Zp,ptBR,ptBZ,ptBphi,ptR0-Rp,ptB0,ptBT0,ptBp0)/Rp;%Safety factor
0075 %
0076 r_lcs = sqrt((ptR(nr,:) - Rp).^2 + (ptZ(nr,:) - Zp).^2);%Minor radius of the last closed magnetic surface
0077 %
0078 smax = 3*(Rp + ap);%maximum chord length = R_hxr+a_hxr
0079 %s = linspace(-smax,+smax,ns)'*ones(1,nc);%for nc (s negative by definition)
0080 s = linspace(-smax,0,ns).'*ones(1,nc);%for nc (s negative by definition)
0081 % Poloidal plot
0082 RR_pol = NaN(nr,nt);
0083 ZZ_pol = NaN(nr,nt);
0084 BR_pol = NaN(nr,nt);
0085 BZ_pol = NaN(nr,nt);
0086 %
0087 for ir = 1:nr_S,
0088     RR_pol(ir,:) = interp1(ppsin,ptR,xpsin_S(ir),method);
0089     ZZ_pol(ir,:) = interp1(ppsin,ptZ,xpsin_S(ir),method);
0090     BR_pol(ir,:) = interp1(ppsin,ptBR,xpsin_S(ir),method);
0091     BZ_pol(ir,:) = interp1(ppsin,ptBZ,xpsin_S(ir),method);
0092 end
0093 %
0094 HHR_pol = linspace(Rp-1.1*r_lcs(round(nt/2)),Rp+1.1*r_lcs(1),1000);
0095 HHZ_pol = Zp + zeros(1,1000);
0096 VVR_pol = Rp + zeros(1,1000);
0097 VVZ_pol = linspace(-1.1*r_lcs(round(nt/4)),1.1*r_lcs(round(3*nt/4)),1000);
0098 % Toroidal plot
0099 XX_tor_min = (Rp-r_lcs(round(nt/2)))*cos(linspace(0,2*pi,1000));
0100 YY_tor_min = (Rp-r_lcs(round(nt/2)))*sin(linspace(0,2*pi,1000));
0101 XX_tor = Rp*cos(linspace(0,2*pi,1000));
0102 YY_tor = Rp*sin(linspace(0,2*pi,1000));
0103 XX_tor_max = (Rp+r_lcs(1))*cos(linspace(0,2*pi,1000));
0104 YY_tor_max = (Rp+r_lcs(1))*sin(linspace(0,2*pi,1000));
0105 VVX = zeros(1,1000);
0106 VVY = linspace(1.1*min(YY_tor_max),1.1*max(YY_tor_max),1000);
0107 HHX = linspace(1.1*min(XX_tor_max),1.1*max(XX_tor_max),1000);
0108 HHY = zeros(1,1000);
0109 %
0110 phi_mfl = linspace(0,2*pi*xq(end),1000);
0111 phi_mfl0 = 0;
0112 theta_mfl_lcs = (phi_mfl - phi_mfl0)/xq(end);%Heliticity relation given by the safety factor
0113 r_mfl_lcs = interp1(theta,r_lcs,theta_mfl_lcs,method);
0114 %
0115 BR_mfl_lcs = interp1(theta,ptBR(end,:),theta_mfl_lcs,method);
0116 Bphi_mfl_lcs = interp1(theta,ptBphi(end,:),theta_mfl_lcs,method);
0117 %
0118 ptBtheta = -ptBR.*sin(ones(nr,1)*theta) + ptBZ.*cos(ones(nr,1)*theta);
0119 Ipdir = sign(ptBtheta);%Btheta < 0 (clockwise poloidal) then Ip is counterclockwise toroidally from the top view (positive trigonometric). May change as function of plasma minor radius (negative shear case)
0120 Ipdir_mfl_lcs = sign(interp1(theta,Ipdir(end,:),theta_mfl_lcs,method));%Sign of Ip on the last closed magnetic surface
0121 %
0122 BX_mfl_lcs = BR_mfl_lcs.*cos(phi_mfl) - Bphi_mfl_lcs.*sin(phi_mfl);
0123 BY_mfl_lcs = BR_mfl_lcs.*sin(phi_mfl) + Bphi_mfl_lcs.*cos(phi_mfl);
0124 %
0125 X_mfl_lcs = (Rp + r_mfl_lcs.*cos(theta_mfl_lcs)).*cos(phi_mfl); 
0126 Y_mfl_lcs = (Rp + r_mfl_lcs.*cos(theta_mfl_lcs)).*sin(phi_mfl); 
0127 Z_mfl_lcs = Zp + r_mfl_lcs.*sin(theta_mfl_lcs);
0128 %
0129 alpha_hxr = ones(ns,1)*hxr.alpha_hxr;%Defined in the indirect frame only (see the documentation)
0130 beta_hxr = ones(ns,1)*hxr.beta_hxr;
0131 R_hxr = ones(ns,1)*hxr.R_hxr;
0132 Z_hxr = ones(ns,1)*hxr.Z_hxr;
0133 %
0134 if any(alpha_hxr(:) <= -pi | alpha_hxr(:) > pi)
0135      error('need -pi < alpha_hxr <= pi')
0136 end
0137 %
0138 if any(beta_hxr(:) < 0 | beta_hxr(:) > pi),
0139     error('need 0 <= beta_hxr <= pi')
0140 end    
0141 %
0142 sZ = s.*cos(beta_hxr);
0143 sh = s.*sin(beta_hxr);
0144 %
0145 sR = sqrt(R_hxr.^2 + sh.^2 + 2*R_hxr.*sh.*cos(alpha_hxr));
0146 sR = sR + eps*(sR == Rp);%to avoid singularities
0147 sZ = Z_hxr + sZ;
0148 %
0149 sZ = sZ + eps*(sZ == Zp);%to avoid singularities
0150 scosphi = (sR.^2 + R_hxr.^2 - sh.^2)./(2*sR.*R_hxr);
0151 scosphi(scosphi < -1) = -1;
0152 scosphi(scosphi > 1) = 1;
0153 %sphi = rem(2*pi - acos(scosphi).*sign(alpha_hxr),2*pi);%Toroidal angle
0154 sphi = -acos(scosphi).*sign(alpha_hxr);%Toroidal angle
0155 spsi = griddata(ptR,ptZ,ppsin(:)*ones(1,nt),sR,sZ);%Flux surface value
0156 %
0157 gX_hxr = sin(beta_hxr).*cos(alpha_hxr);%Chord vector
0158 gY_hxr = sin(beta_hxr).*sin(alpha_hxr);%Chord vector
0159 gZ_hxr = cos(beta_hxr);%Chord vector
0160 %
0161 sr = sqrt((sR - Rp).^2 + (sZ - Zp).^2);%Minor radius along the chords
0162 st = atan((sZ - Zp)./(sR - Rp)) + pi - pi*(sR > Rp).*sign(sZ - Zp);%poloidal angle [0,2*pi]
0163 %
0164 ir_spsi = spsi*0;
0165 for i = nr_S-1:-1:1,
0166     ir_spsi = ir_spsi + (spsi > xpsin_S(i));    
0167 end
0168 %
0169 it_st = st*0;
0170 for i = nt-1:-1:1,
0171     it_st = it_st + (st > theta(i));    
0172 end
0173 %
0174 mask_hxr = zeros(nr_S-1,nt-1);%Location where the distribution function is calculated
0175 for i = 1:nr_S-1,
0176     for j= 1:nt-1,
0177         if ~isempty(find(it_st(:,1:nc) == j & ir_spsi(:,1:nc) == i,1,'first')),
0178             mask_hxr(i,j) = 1;
0179         end
0180     end
0181 end
0182 %
0183 mask = ~isnan(spsi);
0184 for ic = 1:nc,
0185    imask = find(mask(1:end-1,ic)-mask(2:end,ic) == 1);
0186    if length(imask) == 2,%reject the case when chord go through the plasma twice
0187         mask(1:imask(1),ic) = false;
0188    end
0189 end    
0190 %
0191 s(~mask) = NaN;
0192 sR(~mask) = NaN;
0193 sZ(~mask) = NaN; 
0194 sphi(~mask) = NaN; 
0195 sr(~mask) = NaN;
0196 st(~mask) = NaN;
0197 spsi(~mask) = NaN;
0198 ir_spsi(~mask) = NaN;
0199 it_st(~mask) = NaN;
0200 gX_hxr(~mask) = NaN;
0201 gY_hxr(~mask) = NaN;
0202 gZ_hxr(~mask) = NaN;
0203 %
0204 sBR = griddata(ptR,ptZ,ptBR,sR,sZ);
0205 sBZ = griddata(ptR,ptZ,ptBZ,sR,sZ);
0206 sBphi = griddata(ptR,ptZ,ptBphi,sR,sZ);
0207 sB = sqrt(sBR.^2 + sBZ.^2 + sBphi.^2);
0208 sksid = ((sBR.*cos(sphi) - sBphi.*sin(sphi)).*gX_hxr + (-sBR.*sin(sphi) + sBphi.*cos(sphi)).*gY_hxr + sBZ.*gZ_hxr)./sB;%cosine between magnetic field line and chord direction (towards the detector)
0209 %
0210 vschords = zeros(16,ns,nc);
0211 vschords(1,:,:) = s;
0212 vschords(2,:,:) = sR;
0213 vschords(3,:,:) = sZ;
0214 vschords(4,:,:) = sphi;
0215 vschords(5,:,:) = sr;
0216 vschords(6,:,:) = st;
0217 vschords(7,:,:) = spsi;
0218 vschords(8,:,:) = ir_spsi;
0219 vschords(9,:,:) = it_st;
0220 vschords(10,:,:) = sksid;
0221 vschords(11,:,:) = sBR;
0222 vschords(12,:,:) = sBZ;
0223 vschords(13,:,:) = sBphi;
0224 vschords(14,:,:) = gX_hxr;
0225 vschords(15,:,:) = gY_hxr;
0226 vschords(16,:,:) = gZ_hxr;
0227 %
0228 % Structure for display
0229 %
0230 Zbremchord.ptrho = ptrho;
0231 Zbremchord.xq = xq;
0232 Zbremchord.sR = sR;
0233 Zbremchord.sZ = sZ;
0234 Zbremchord.RR_pol = RR_pol;
0235 Zbremchord.ZZ_pol = ZZ_pol;
0236 Zbremchord.HHR_pol = HHR_pol;
0237 Zbremchord.HHZ_pol = HHZ_pol;
0238 Zbremchord.VVR_pol = VVR_pol;
0239 Zbremchord.VVZ_pol = VVZ_pol;
0240 Zbremchord.ptR = ptR;
0241 Zbremchord.ptZ = ptZ;
0242 Zbremchord.BR_pol = BR_pol;
0243 Zbremchord.BZ_pol = BZ_pol;
0244 Zbremchord.X_mfl_lcs = X_mfl_lcs;
0245 Zbremchord.Y_mfl_lcs = Y_mfl_lcs;
0246 Zbremchord.Z_mfl_lcs = Z_mfl_lcs;
0247 Zbremchord.BX_mfl_lcs = BX_mfl_lcs;
0248 Zbremchord.BY_mfl_lcs = BY_mfl_lcs;
0249 Zbremchord.Ipdir_mfl_lcs = Ipdir_mfl_lcs;
0250 Zbremchord.sR = sR;
0251 Zbremchord.sphi = sphi;
0252 Zbremchord.XX_tor_min = XX_tor_min;
0253 Zbremchord.YY_tor_min = YY_tor_min;
0254 Zbremchord.XX_tor = XX_tor;
0255 Zbremchord.YY_tor = YY_tor;
0256 Zbremchord.XX_tor_max = XX_tor_max;
0257 Zbremchord.YY_tor_max = YY_tor_max;
0258 Zbremchord.VVX = VVX;
0259 Zbremchord.VVY = VVY;
0260 Zbremchord.HHX = HHX;
0261 Zbremchord.HHY = HHY;
0262 Zbremchord.xpsin_S = xpsin_S;
0263 Zbremchord.vschords = vschords;
0264 %