Calculation of the critical poloidal angle value theta_star for runaway avalanches INPUT: - equil: magnetic equilibrium structure - equil_fit: structure of the vectorial for of the magnetic equilibrium - gridDKE: grid structure - Zmomcoef: momentum dynamics structure - radialDKE: radial grid structure OUTPUT: - XXthetastar1: poloidal angle value theta_star between 0 and pi [np,nmhu,nr] - XXthetastar2: poloidal angle value theta_star between pi and 2*pi [np,nmhu,nr] by Y. PEYSSON (CEA/DRFC, yves.peysson@cea.fr) and D. SCHMITT (CEA/DRFC, damien.schmitt@cea.fr)
0001 function [XXthetastar1,XXthetastar2] = thetastar_dke_ds(equil,equil_fit,gridDKE,Zmomcoef,radialDKE) 0002 % 0003 % 0004 % Calculation of the critical poloidal angle value theta_star for runaway avalanches 0005 % 0006 % INPUT: 0007 % 0008 % - equil: magnetic equilibrium structure 0009 % - equil_fit: structure of the vectorial for of the magnetic equilibrium 0010 % - gridDKE: grid structure 0011 % - Zmomcoef: momentum dynamics structure 0012 % - radialDKE: radial grid structure 0013 % 0014 % OUTPUT: 0015 % 0016 % - XXthetastar1: poloidal angle value theta_star between 0 and pi [np,nmhu,nr] 0017 % - XXthetastar2: poloidal angle value theta_star between pi and 2*pi [np,nmhu,nr] 0018 % 0019 % by Y. PEYSSON (CEA/DRFC, yves.peysson@cea.fr) and D. SCHMITT (CEA/DRFC, damien.schmitt@cea.fr) 0020 % 0021 [equil_fit] = fitequil_yp(equil);%Vectorial form of the magnetic equilibrium 0022 % 0023 npn = size(gridDKE.XXmhu,1); 0024 nmhu = size(gridDKE.XXmhu,2); 0025 nr = size(gridDKE.XXmhu,3); 0026 ntheta = length(equil.theta); 0027 % 0028 xpsin = radialDKE.xpsin_f; 0029 % 0030 XXB0 = permute(repmat(sqrt(equil.ptBx(:,1).^2 + equil.ptBy(:,1).^2 + equil.ptBPHI(:,1).^2),[1,npn,nmhu]),[2,3,1]); 0031 XXBmax = permute(repmat(sqrt(equil.ptBx(:,(ntheta+1)/2).^2 + equil.ptBy(:,(ntheta+1)/2).^2 + equil.ptBPHI(:,(ntheta+1)/2).^2),[1,npn,nmhu]),[2,3,1]); 0032 XXBstar = 2*XXB0./(1 + Zmomcoef.XXgamma.^2)./(1 - gridDKE.XXmhu.^2); 0033 % 0034 % *************** Put criterion that XXBstar >= XXB0 and XXBstar <= XXBmax********************** 0035 % 0036 options = optimset('Display','off'); 0037 % 0038 for ip = 1:npn, 0039 for im = 1:nmhu, 0040 for ir = 1:nr, 0041 [thetas1,Bs1,flags1] = fsolve(@(theta) Bfit_dke_ds(theta,xpsin(ir),equil,equil_out,XXBstar(ip,im,ir)),pi/2);%Solution between [0,pi] 0042 [thetas2,Bs2,flags2] = fsolve(@(theta) Bfit_dke_ds(theta,xpsin(ir),equil,equil_out,XXBstar(ip,im,ir)),3*pi/2);%Solution between [pi,2*pi] 0043 if flags1 == 1, XXthetastar1(ip,im,ir) = thetas1; else XXthetastar1(ip,im,ir) = NaN; end 0044 if flags2 == 1, XXthetastar2(ip,im,ir) = thetas2; else XXthetastar2(ip,im,ir) = NaN; end 0045 end 0046 end 0047 end 0048 % 0049 0050 0051