rz2pt_jd

 This function inverts the magnetic geometry data.

 by J. Decker (jodecker@alum.mit.edu) 11/02/04

0001 function [ptx,pty,ptdpsindx,ptdpsindy,psin,theta] = rz2pt_jd(x,y,xypsin,psin,theta,method,display_mode,p_opt,filename,opt)
0002 %
0003 % This function inverts the magnetic geometry data.
0004 %
0005 % by J. Decker (jodecker@alum.mit.edu) 11/02/04
0006 %
0007 if nargin < 10
0008     opt.sfac = 10;% factor for (x,y) supergrid
0009     opt.axenforce = true;% option to enforce center position and psi value in supergrid
0010 end
0011 if nargin < 9
0012     filename = 'tokamak';
0013 end
0014 if nargin < 8
0015     p_opt = -1;
0016 end
0017 if nargin < 7
0018     display_mode = 0;
0019 end
0020 if nargin < 6
0021     method = 'linear';
0022 end
0023 if nargin < 5
0024     theta = 65;
0025 end
0026 if nargin < 4
0027     psin = 101;
0028 end
0029 %
0030 nt = length(theta);
0031 if nt == 1
0032     nt = theta;
0033     theta = linspace(0,2*pi,nt);
0034 end
0035 %
0036 np = length(psin);
0037 if np == 1
0038     np = psin;
0039     psin = linspace(0,1,np);
0040 end
0041 %
0042 psin_lcfs = 0.999;%0.998;%
0043 %
0044 nx=  length(x);
0045 ny = length(y);
0046 %
0047 % interpolation for better griddata
0048 %
0049 xmin = min(x);
0050 xmax = max(x);
0051 ymin = min(y);
0052 ymax = max(y);
0053 %
0054 snx = opt.sfac*(nx-1) + 1;
0055 sny = opt.sfac*(ny-1) + 1;
0056 %
0057 sx = linspace(xmin,xmax,snx);
0058 sy = linspace(ymin,ymax,sny);
0059 %
0060 % even finer grid near the center
0061 %
0062 %ix = find(abs(sx) == min(abs(sx)));
0063 %iy = find(abs(sy) == min(abs(sy)));
0064 %
0065 %sx = sort([sx([1:ix-2,ix+2:end]),linspace(sx(ix - 1),sx(ix + 1),2*opt.sfac+1)]);
0066 %sy = sort([sy([1:iy-2,iy+2:end]),linspace(sy(iy - 1),sy(iy + 1),2*opt.sfac+1)]);
0067 %
0068 %
0069 %snx = length(sx)
0070 %sny = length(sy)
0071 %
0072 if opt.axenforce,
0073     ix = find(abs(sx) == min(abs(sx)));
0074     iy = find(abs(sy) == min(abs(sy)));
0075     %
0076     if abs(sx(ix)) > 10*eps || abs(sy(iy)) > 10*eps,
0077         disp(' ---> WARNING: location of magnetic axis enforced in rz2pt_jd.m')
0078     end
0079     sx(ix) = 0;
0080     sy(iy) = 0;
0081 end
0082 %
0083 sxypsin = interp2(x,y,xypsin.',sx.',sy,method).';
0084 %
0085 if opt.axenforce,
0086     if abs(sxypsin(ix,iy)) > 10*eps,
0087         disp(' ---> WARNING: value of psi at magnetic axis enforced in rz2pt_jd.m')
0088     end
0089     sxypsin(ix,iy) = 0;
0090     %
0091 end
0092 %
0093 % calculation of grad psi
0094 %
0095 hx = (sx(1:snx-1) + sx(2:snx))/2;
0096 hy = (sy(1:sny-1) + sy(2:sny))/2;
0097 %
0098 hxyx = repmat(hx.',[1,sny]);
0099 hxyy = repmat(hy,[snx,1]);
0100 %
0101 hxypsin_x = interp1(sx,sxypsin,hx,method);
0102 hxypsin_y = interp1(sy,sxypsin.',hy,method).';
0103 %
0104 sxydpsindx = zeros(snx,sny);
0105 sxydpsindy = zeros(snx,sny);
0106 %
0107 sxydpsindx(2:snx-1,:) = diff(hxypsin_x)./diff(hxyx);
0108 sxydpsindy(:,2:sny-1) = diff(hxypsin_y.').'./diff(hxyy.').';
0109 %
0110 sxydpsindx(1,:) = (hxypsin_x(1,:) - sxypsin(1,:))/(hx(1) - sx(1));
0111 sxydpsindy(:,1) = (hxypsin_y(:,1) - sxypsin(:,1))/(hy(1) - sy(1));
0112 %
0113 sxydpsindx(snx,:) = (sxypsin(snx,:) - hxypsin_x(snx-1,:))/(sx(snx) - hx(snx-1));
0114 sxydpsindy(:,sny) = (sxypsin(:,sny) - hxypsin_y(:,sny-1))/(sy(sny) - hy(sny-1));
0115 %
0116 if opt.axenforce,
0117     sxydpsindx(ix,iy) = 0;
0118     sxydpsindy(ix,iy) = 0;
0119 end
0120 %
0121 % change to (r,theta) grid
0122 %
0123 sxyx = repmat(sx.',[1,sny]);
0124 sxyy = repmat(sy,[snx,1]);
0125 %
0126 sxyr = sqrt(sxyx.^2 + sxyy.^2);
0127 sxyx(sxyx == 0) = - eps;
0128 sxytheta = atan(sxyy./sxyx) + pi*(sxyx < 0) + 2*pi*(sxyx > 0 & sxyy < 0);
0129 %
0130 % last closed flux-surface
0131 %
0132 [xc,yc] = contlimit_jd(sx,sy,sxypsin',[-1,psin_lcfs]);
0133 %
0134 % selects the plasma among various closed systems
0135 %
0136 for ic = 1:length(xc),
0137     nc(ic) = length(xc{ic});
0138 end
0139 ic = find(nc == max(nc)); 
0140 sxa = xc{ic};
0141 sya = yc{ic};
0142 %
0143 % change to (r,theta) grid
0144 %
0145 sra = sqrt(sxa.^2 + sya.^2);
0146 sxa(sxa == 0) = - eps;
0147 sthetaa = atan(sya./sxa) + pi*(sxa < 0) + 2*pi*(sxa > 0 & sya < 0);
0148 [sthetaa,ia] = sort(sthetaa);
0149 sra = sra(ia);
0150 maska = [1,find(diff(sthetaa) ~= 0)+1];
0151 sthetaa = sthetaa(maska);
0152 sra = sra(maska);
0153 sna = length(sthetaa);
0154 sthetaa = [sthetaa(sna) - 2*pi,sthetaa,sthetaa(1) + 2*pi];
0155 sra = [sra(sna),sra,sra(1)];
0156 %
0157 % redefine free space values for better interpolations
0158 %
0159 sXtheta = reshape(sxytheta,[snx*sny,1]);
0160 sXr = reshape(sxyr,[snx*sny,1]);
0161 sXpsin = reshape(sxypsin,[snx*sny,1]);
0162 sXdpsindx = reshape(sxydpsindx,[snx*sny,1]);
0163 sXdpsindy = reshape(sxydpsindy,[snx*sny,1]);
0164 %
0165 sXrmax = interp1(sthetaa,sra,sXtheta,method);
0166 sXmask_in = sXr <= sXrmax;
0167 %
0168 sXtheta = sXtheta(sXmask_in);
0169 sXr = sXr(sXmask_in);
0170 sXpsin = sXpsin(sXmask_in);
0171 sXdpsindx = sXdpsindx(sXmask_in);
0172 sXdpsindy = sXdpsindy(sXmask_in);
0173 %
0174 % Remove central point (which is undefined)
0175 %
0176 smask0 = (sXr ~= 0);
0177 %
0178 sXtheta = sXtheta(smask0);
0179 sXr = sXr(smask0);
0180 sXpsin = sXpsin(smask0);
0181 sXdpsindx = sXdpsindx(smask0);
0182 sXdpsindy = sXdpsindy(smask0);
0183 %
0184 % For better theta calculations
0185 %
0186 %sXtheta = [sXtheta;ttheta'];
0187 %sXr = [sXr;zeros(nt,1)];
0188 %sXpsin = [sXpsin;zeros(nt,1)];
0189 %
0190 smaskp = find(sXtheta >= 0 & sXtheta < pi/8);
0191 sXtheta = [sXtheta;2*pi+sXtheta(smaskp)];
0192 sXr = [sXr;sXr(smaskp)];
0193 sXpsin = [sXpsin;sXpsin(smaskp)];
0194 sXdpsindx = [sXdpsindx;sXdpsindx(smaskp)];
0195 sXdpsindy = [sXdpsindy;sXdpsindy(smaskp)];
0196 %
0197 smaskm = find(sXtheta < 2*pi & sXtheta > 2*pi - pi/8);
0198 sXtheta = [sXtheta;sXtheta(smaskm) - 2*pi];
0199 sXr = [sXr;sXr(smaskm)];
0200 sXpsin = [sXpsin;sXpsin(smaskm)];
0201 sXdpsindx = [sXdpsindx;sXdpsindx(smaskm)];
0202 sXdpsindy = [sXdpsindy;sXdpsindy(smaskm)];
0203 %
0204 hpsin = psin(2:np-1);
0205 %
0206 gdmethod = 'cubic';
0207 %
0208 warning off
0209 disp(' ---> WARNING: Duplicate x-y data points detected: using average values for duplicate points in rz2pt_jd.m')
0210 hptr = griddata(sXpsin,sXtheta,sXr,hpsin.',theta,gdmethod).';
0211 tra = interp1(sthetaa,sra,theta,method);
0212 ptr = [zeros(1,nt);hptr;tra];
0213 %
0214 pttheta = repmat(theta,[np,1]);
0215 %
0216 ptx = ptr.*cos(pttheta);
0217 pty = ptr.*sin(pttheta);
0218 %
0219 %sXR = reshape(srzR,[snr*snz,1]);
0220 %sXZ = reshape(srzZ,[snr*snz,1]);
0221 %sXR = sXR(sXmask_in);
0222 %sXZ = sXZ(sXmask_in);
0223 %sXR = sXR(~smask0);
0224 %sXZ = sXZ(~smask0);
0225 %sXR = [sXR;sXR(smaskp)];
0226 %sXZ = [sXZ;sXZ(smaskp)];
0227 %sXR = [sXR;sXR(smaskm)];
0228 %sXZ = [sXZ;sXZ(smaskm)];
0229 %p2tR = griddata(sXpsin,sXtheta,sXR,p2psin',ttheta,gdmethod)';
0230 %p2tZ = griddata(sXpsin,sXtheta,sXZ,p2psin',ttheta,gdmethod)';
0231 %
0232 %h2psin = [hpsin(1:end-1),linspace(hpsin(end),1,opt.sfac)];
0233 h2psin = hpsin;
0234 %
0235 disp(' ---> WARNING: Duplicate x-y data points detected: using average values for duplicate points in rz2pt_jd.m')
0236 h2ptdpsindx = griddata(sXpsin,sXtheta,sXdpsindx,h2psin.',theta,gdmethod).';
0237 disp(' ---> WARNING: Duplicate x-y data points detected: using average values for duplicate points in rz2pt_jd.m')
0238 h2ptdpsindy = griddata(sXpsin,sXtheta,sXdpsindy,h2psin.',theta,gdmethod).';
0239 warning on
0240 %
0241 %p2tdpsindr
0242 %
0243 %hXR = reshape(hrzR,[hnr*hnz,1]);
0244 %hXZ = reshape(hrzZ,[hnr*hnz,1]);
0245 %hXR = hXR(hXmask_in);
0246 %hXZ = hXZ(hXmask_in);
0247 %hXR = [hXR;hXR(hmaskp)];
0248 %hXZ = [hXZ;hXZ(hmaskp)];
0249 %hXR = [hXR;hXR(hmaskm)];
0250 %hXZ = [hXZ;hXZ(hmaskm)];
0251 %p2tR = griddata(hXpsin,hXtheta,hXR,p2psin',ttheta,gdmethod)';
0252 %p2tZ = griddata(hXpsin,hXtheta,hXZ,p2psin',ttheta,gdmethod)';
0253 %
0254 ipmax = max(find(sum(isnan(h2ptdpsindx.') | isnan(h2ptdpsindy.')) == 0));
0255 mask = 1:ipmax;
0256 %
0257 interp1method = 'pchip';%'cubic' method not recognized in interp1
0258 %
0259 ptdpsindx = [zeros(1,nt);interp1(h2psin(mask),h2ptdpsindx(mask,:),[hpsin,1],interp1method,'extrap')];
0260 ptdpsindy = [zeros(1,nt);interp1(h2psin(mask),h2ptdpsindy(mask,:),[hpsin,1],interp1method,'extrap')];
0261 %
0262 if display_mode >= 2,
0263     %
0264     style = '-';
0265     marker = 'none';
0266     color1 = 'r';
0267     color2 = 'b';
0268     width1 = 2;
0269     width2 = 0.5;
0270     siz = 20;
0271     red = 0.9;
0272     %
0273     ap = ptx(np,1);%Plasma minor radius on the LFS midplane
0274     prho = ptx(:,1).'/ap;
0275     %
0276     srho = linspace(0,1,21);
0277     stheta = linspace(0,2*pi,33);
0278     sx = interp1(prho,ptx,srho,method);  
0279     sy = interp1(prho,pty,srho,method); 
0280     stx = interp1(theta,ptx',stheta,method)';  
0281     sty = interp1(theta,pty',stheta,method)';    
0282     %
0283     xmin = min(min(ptx));
0284     xmax = max(max(ptx));
0285     ymin = min(min(pty));
0286     ymax = max(max(pty));
0287     %
0288     figure_yp('','Poloidal view - MHD equilibrium','clf');
0289     %
0290     [ax] = graph1D_jd(sx',sy',0,0,'','','',NaN,[xmin,xmax],[ymin,ymax],style,marker,color1,width1,siz,gca,red);pause(0.5)
0291     [ax] = graph1D_jd(stx,sty,0,0,'x','y','',NaN,[xmin,xmax],[ymin,ymax],style,marker,color2,width2);pause(0.5)
0292     axis('equal');  
0293     axis([xmin,xmax,ymin,ymax]);
0294     print_jd(p_opt,[filename,'_rhotheta']);
0295     %
0296     figure_yp('','Poloidal view - MHD equilibrium','clf');
0297     %
0298     sspsin = linspace(0,1,21);
0299     sx = interp1(sqrt(psin),ptx,sspsin,method);  
0300     sy = interp1(sqrt(psin),pty,sspsin,method); 
0301     %
0302     [ax] = graph1D_jd(sx',sy',0,0,'x','y','',NaN,[xmin,xmax],[ymin,ymax],style,marker,color1,width1,siz,gca,red);pause(0.5)
0303     axis('equal');  
0304     axis([xmin,xmax,ymin,ymax]);
0305     print_jd(p_opt,[filename,'_spsin']);
0306     %
0307 end