thetaconemj3

 calculation of ripple losscone angles

0001 % calculation of ripple losscone angles
0002 %
0003 Ip   = Ip*1000000;
0004 alfap1 = 4;
0005 B0R0   = B0*R0;
0006 CBP    = 2e-7*Ip;
0007 NBOB   = 18;
0008 B0     = B0R0/R0;
0009 Z0     = 0;
0010 a2     = a*a;
0011 Raxe   = R0 + D0;
0012 c0     = D0/a2;
0013 rc2    = c0*c0;
0014 dc2    = 2*rc2;
0015 qc2    = 4*rc2;
0016 %
0017 ag         = 1.2;
0018 RR         = 1.54:0.025:1.54+ag*11/8;
0019 ZZ         = 0:0.025:1.2;
0020 [R,Z]      = meshgrid(RR,ZZ);
0021 zerR       = zeros(size(R));
0022 X          = (R-Raxe);
0023 Y          = 1-2*c0*X;
0024 rs2        = (Y-sqrt(Y.*Y-qc2*(X.^2+(Z-Z0).^2)))/dc2;
0025 rs         = sqrt(rs2);
0026 cr2        = max(zerR,1-rs2/a2);
0027 BpM        = CBP*(1-cr2.^alfap1)./rs;
0028 D          = D0*(1-rs2/a2);
0029 Dprime     = -2*D0*rs/a2;
0030 cteta      = (R-R0-D)./rs;
0031 Bp         = BpM.*(R0+D)./(R.*(1+Dprime.*cteta));
0032 %
0033 BpR        = Bp.*((Z-Z0)./rs);
0034 BpZ        = -Bp.*cteta;
0035 BF         = B0R0./R;
0036 %
0037 teta       = (0:0.005:1)*2*pi;
0038 Rbord      = R0+a*cos(teta);
0039 Zbord      = Z0+a*sin(teta);
0040 %
0041 ripple        = deltacalmj3(R,Z);
0042 alfastar      = abs(BpR./(NBOB*ripple.*BF));
0043 
0044 rhos  = rho*a;
0045 Rcs   = R0+D0*(1-(rhos/a).^2);
0046 Zcs   = zeros(size(rhos));
0047 Rmini = Rcs-rhos;
0048 Rmaxi = Rcs+rhos;
0049 
0050 cc=contour(RR,ZZ,alfastar,[-1 1]);clf;close
0051 Rastar   = cc(1,2:cc(2,1));
0052 Zastar   = cc(2,2:cc(2,1));
0053 Rastar(end+1) = Rastar(1);
0054 Zastar(end+1) = Zastar(1);
0055 centreR=mean(Rastar);centreZ=mean(Zastar);
0056 angastar=unwrap(angle((Rastar-centreR)+i*(Zastar-centreZ)));
0057 
0058 newang      = linspace(angastar(1),angastar(end),1000);
0059 newang(end) = [];
0060 nRastar     = interp1(angastar,Rastar,newang,'linear');
0061 nZastar     = interp1(angastar,Zastar,newang,'linear');
0062 nRastar(isnan(nRastar)) = -ones(size(nRastar(isnan(nRastar))));
0063 nZastar(isnan(nZastar)) = -ones(size(nZastar(isnan(nZastar))));
0064 
0065 R           = linspace(Rmaxi-1e-6,Rmini+1e-6,1000);
0066 Z           = sqrt(rhos^2-(R-Rcs).^2);
0067 nZ          = sqrt(rhos^2-(nRastar-Rcs).^2);
0068 ind         = find(~imag(nZ));
0069 if isempty(ind)
0070   disp(['WARNING: no intersection between the good confinement area and the flux surface at rho=',num2str(rhos,4)])
0071   mulim           = zeros(size(rhos));
0072   thetabanane     = zeros(size(rhos));
0073   thetaloss       = zeros(size(rhos));
0074   thetaZBC        = zeros(size(rhos));
0075   return
0076 
0077 end
0078 dis         = nZ-nZastar;
0079 dis1        = cat(2,dis(ind),0);
0080 dis2        = cat(2,0,dis(ind));
0081 indval         = find((dis1.*dis2) <=0);
0082 %if length(indval) > 1,
0083 %   indval(1)      = [];
0084 %   indval(end)    = [];
0085 %else
0086 %   indval =[];
0087 %end
0088 indval(1)      = [];
0089 indval(end)    = [];
0090 indRb        = ind(indval-1);
0091 Rb           = nRastar(indRb);
0092 Zb           = nZastar(indRb);
0093 [Rb,iRb]     = max(Rb);
0094 Zb           = Zb(iRb);
0095 
0096 if rho>0.7 & length(indRb) ==0 
0097    Rb = 10;
0098    Zb = 10;
0099 end
0100 
0101 mulim           = zeros(size(rhos));
0102 thetabanane     = zeros(size(rhos));
0103 thetaloss       = zeros(size(rhos));
0104 thetaZBC        = zeros(size(rhos));
0105 
0106 mubanane        = sqrt(2*rhos./(R0+rhos));
0107 thetabanane     = asin(mubanane)*180/pi;
0108 
0109 ind1            = find(Rb>0&Rb<10&(Rb-R0)>0);
0110 thetaZBC(ind1)  = atan(Zb(ind1)./(Rb(ind1)-R0));
0111 ind2            = find(Rb>0&Rb<10&(Rb-R0)<=0);
0112 thetaZBC(ind2)  = 0.5*pi+atan((R0-Rb(ind2))./Zb(ind2));
0113 ind3            = find(Rb==10);
0114 thetaZBC(ind3)  = pi*ones(1,length(ind3));
0115 thetaZBC        = thetaZBC*180/pi;
0116 
0117 ind             = [ind1 ind2];
0118 mulim(ind)      = sqrt(rhos(ind).*(1-cos(thetaZBC(ind)*pi/180))./(rhos(ind)+R0));
0119 mulim(ind3)     = sqrt(2*rhos(ind3)./(R0+rhos(ind3)));
0120 thetaloss(ind)  = asin(mulim(ind))*180/pi;
0121 thetaloss(ind3) = asin(mulim(ind3))*180/pi;