equil_magnetic_AUG_jd

0001 function [equil,prho,pspsin] = equil_magnetic_AUG_jd(external,opt)
0002 %
0003 if nargin < 2 || ~isstruct(opt),
0004     opt = struct;
0005 end
0006 %
0007 if ~isfield(opt,'nlev'),
0008     opt.nlev = 100;
0009 end
0010 if ~isfield(opt,'axis'),
0011     opt.axis = true;%recalculatethe axis position
0012 end
0013 if ~isfield(opt,'np'),
0014     opt.np = 100;
0015 end
0016 if ~isfield(opt,'nt'),
0017     opt.nt = 65;
0018 end
0019 %
0020 %
0021 R = external.R;
0022 Z = external.Z;
0023 BR = external.BR;
0024 BZ = external.BZ;
0025 Bphi = external.Bphi;
0026 psi = external.psi;
0027 %
0028 if all(isfield(external,{'Rp','Zp','psi0','psia'})),
0029     Rp = external.Rp;
0030     Zp = external.Zp;
0031     psip = external.psi0;
0032     psia = external.psia;
0033 else
0034     %
0035     [Rc,Zc,~,lev] = contlimit_jd(R,Z,psi.',opt.nlev);
0036     %
0037     nc = length(lev);
0038     lc = zeros(1,nc);
0039     le = zeros(1,nc);
0040     for ic = 1:nc,
0041         lc(ic) = sum(sqrt((Rc{ic}(2:end) - Rc{ic}(1:end-1)).^2 + (Zc{ic}(2:end) - Zc{ic}(1:end-1)).^2));
0042         la(ic) = (max(Rc{ic}) - min(Rc{ic})).*(max(Zc{ic}) - min(Zc{ic}));
0043         le(ic) = sqrt((Rc{ic}(1) - Rc{ic}(end)).^2 + (Zc{ic}(1) - Zc{ic}(end)).^2);
0044     end
0045     %
0046     icmask = find(lc > 0 & la > 0 & le == 0);
0047     %
0048     % figure
0049     % for iic = 1:length(icmask),
0050     %     plot(Rc{icmask(iic)},Zc{icmask(iic)}),hold on
0051     % end
0052     % axis('equal')
0053     ic_lcfs = icmask(lc(icmask) == max(lc(icmask)));
0054     psia = lev(ic_lcfs);
0055     %
0056     Rmin_lcfs = min(Rc{ic_lcfs});
0057     Rmax_lcfs = max(Rc{ic_lcfs});
0058     Zmin_lcfs = min(Zc{ic_lcfs});
0059     Zmax_lcfs = max(Zc{ic_lcfs});
0060     %
0061     % reject contours outside of lcf
0062     for iic = length(icmask):-1:1,
0063         Rav = mean(Rc{icmask(iic)});
0064         Zav = mean(Zc{icmask(iic)});
0065         if Rav < Rmin_lcfs || Rav > Rmax_lcfs || Zav < Zmin_lcfs || Zav > Zmax_lcfs,
0066             icmask(iic) = [];
0067         end
0068     end
0069     %
0070     ic_ifs = icmask(lc(icmask) == min(lc(icmask)));
0071     psi_ifs = lev(ic_ifs);
0072     %
0073     Rmin_ifs = min(Rc{ic_ifs});
0074     Rmax_ifs = max(Rc{ic_ifs});
0075     Zmin_ifs = min(Zc{ic_ifs});
0076     Zmax_ifs = max(Zc{ic_ifs});
0077     %
0078     icext = find((lc == 0 | la == 0) & le == 0);
0079     %
0080     next = length(icext);
0081     mask_ext = false(1,next);
0082     mask_lcfs = false(1,next);
0083     %
0084     for iic = 1:next,
0085         Rcext = mean(Rc{icext(iic)});
0086         Zcext = mean(Zc{icext(iic)});
0087         mask_ext(iic) = Rcext < Rmax_ifs && Rcext > Rmin_ifs && Zcext < Zmax_ifs && Zcext > Zmin_ifs;
0088         mask_lcfs(iic) = Rcext < Rmax_lcfs && Rcext > Rmin_lcfs && Zcext < Zmax_lcfs && Zcext > Zmin_lcfs;
0089     end
0090     %
0091     if sum(mask_lcfs) == 0,%no point represent plasma center
0092         %
0093         psip_orig = psi_ifs;
0094         Rp_orig = (Rmin_ifs + Rmax_ifs)/2;
0095         Zp_orig = (Zmin_ifs + Zmax_ifs)/2;
0096         %
0097     else
0098         %
0099         if sum(mask_ext) == 1,
0100             %
0101             psip_orig = lev(icext(mask_ext));
0102             Rp_orig =  mean(Rc{icext(mask_ext)});
0103             Zp_orig =  mean(Zc{icext(mask_ext)});
0104             %
0105         else
0106             error('cannot identify plasma major radius')
0107             %
0108         end
0109     end
0110     if opt.axis,
0111         sR = linspace(Rmin_ifs,Rmax_ifs,101);
0112         sZ = linspace(Zmin_ifs,Zmax_ifs,111);
0113         %
0114         spsi = interp2(R,Z,psi.',sR.',sZ,'spline').';
0115         %
0116         sigma = sign(psia - psi_ifs);
0117         %
0118         [isR,isZ] = find(sigma*spsi == min(min(sigma*spsi)));
0119         Rp = sR(isR);
0120         Zp = sZ(isZ);
0121         psip = spsi(isR,isZ);
0122         %
0123         disp(['Old value of psip: ',num2str(psip_orig),', new value: ',num2str(psip)])
0124         disp(['Old value of Rp: ',num2str(Rp_orig),', new value: ',num2str(Rp)])
0125         disp(['Old value of Zp: ',num2str(Zp_orig),', new value: ',num2str(Zp)])
0126     else
0127         Rp = Rp_orig;
0128         Zp = Zp_orig;
0129         psip = psip_orig;
0130     end
0131 end
0132 %
0133 x = R - Rp;
0134 y = Z - Zp;
0135 %
0136 psi = psi - psip;
0137 psia = psia - psip;
0138 xypsin = psi/psia;
0139 %
0140 [ptx,pty,~,~,psin,theta] = rz2pt_jd(x,y,xypsin,opt.np,opt.nt,'spline',0);
0141 %
0142 ap = ptx(opt.np,1);%Plasma minor radius on the LFS midplane
0143 psi_apRp = psin*psia*ap/Rp;%The poloidal flux coordinate is normalized by the aspect ratio
0144 %
0145 ptBx = interp2(x,y,BR.',ptx,pty,'spline');
0146 ptBy = interp2(x,y,BZ.',ptx,pty,'spline');
0147 ptBPHI = interp2(x,y,Bphi.',ptx,pty,'spline');
0148 %
0149 prho = ptx(:,1).'/ap;
0150 pspsin = sqrt(psin);
0151 %
0152 equil.psi_apRp = psi_apRp;
0153 equil.theta = theta;
0154 equil.ptx = ptx;
0155 equil.pty = pty;
0156 equil.ptBx = ptBx;
0157 equil.ptBy = ptBy;
0158 equil.ptBPHI = ptBPHI;
0159 equil.Rp = Rp;
0160 equil.Zp = Zp;
0161 equil.equi_src = external;
0162 %
0163 
0164