Description of interp_theta

0001 function [Ytheta1,Ytheta2,YF1,YF2] = interp_theta_jd(theta,XBx,XBy,XBz,YBb,XF)
0002 %
0003 % This function calculates the positions Ytheta1 and Ytheta2 located respectively between two grid
0004 % positions Xtheta(i1a) and Xtheta(i1b) for Ytheta1, and Xtheta(i2a) and Xtheta(i2b) for Ytheta2,
0005 % obtained by interpolation with the condition B(theta) = YBb, meaning such that, for j=1,2 :
0006 %
0007 % - XBx(ija)*alpha + XBx(ijb)*(1-alpha) := Bxj
0008 % - XBy(ija)*alpha + XBy(ijb)*(1-alpha) := Byj
0009 % - XBz(ija)*alpha + XBz(ijb)*(1-alpha) := Bzj
0010 % - sqrt(Bxj^2 + Byj^2 + Bzj^2) := Bj = YBb
0011 %
0012 % This function is restricted to the case where XB=sqrt(XBx^2 + XBy^2 + XBz^2) has only one maximum
0013 % and one minimum, such that the equation B(theta) = YBb has either 2 solutions or none.
0014 % The function also interpolates the data XF at the theta location
0015 %
0016 %   INPUT :
0017 %
0018 %       - theta : poloidal grid [1,nt]. Must be periodic with Xtheta(nt) = Xtheta(1)
0019 %       - XBx : x component of the magnetic field [nt,np]. Must be periodic with XBx(nt,:) = XBx(1,:)
0020 %       - XBy : y component of the magnetic field [nt,np]. Must be periodic with XBy(nt,:) = XBy(1,:)
0021 %       - XBz : z component of the magnetic field [nt,np]. Must be periodic with XBz(nt,:) = XBz(1,:)
0022 %       - YBb : magnetic field amplitude at position theta [np,nb].
0023 %       - XF : data to be interpolated (optional) [nt,np,nf,...].
0024 %
0025 %   OUTPUT :
0026 %
0027 %       - Ytheta1 : interpolated position, solution 1 [np,nb]. Ytheta1[ip,iY] = NaN if there no solution.
0028 %       - Ytheta2 : interpolated position, solution 2 [np,nb]. Ytheta2[ip,iY] = NaN if there no solution.
0029 %       - XF : data interpolated at theta = theta1 (optional) [np,nb,nf,...].
0030 %       - XF : data interpolated at theta = theta2 (optional) [np,nb,nf,...].
0031 %
0032 % by J. Decker (joan.decker@cea.fr) and Y. Peysson (yves.peysson@cea.fr)
0033 %
0034 %
0035 if nargin < 5,
0036     error('Not enough input arguments');
0037 end
0038 %
0039 theta = theta(:);
0040 nt = length(theta);
0041 if nt < 3,
0042     error('The number of theta points must be at least 3');
0043 end
0044 %
0045 sBx = size(XBx);
0046 sBy = size(XBy);
0047 sBz = size(XBz);
0048 sBb = size(YBb);
0049 %
0050 if length(sBx) ~= 2 | sBx(1) ~= nt,
0051     error('Input variables dimensions mismatch');
0052 else
0053     np = sBx(2);
0054 end
0055 %
0056 if length(sBy) ~= 2 | length(sBz) ~= 2 | sBy ~= sBx | sBz ~= sBx | sBb(1) ~= np,
0057     error('Input variables dimensions mismatch');
0058 end
0059 %
0060 nb = size(YBb,2);
0061 %
0062 if nargin < 6,
0063     XF = zeros(np,nb,0);
0064 end
0065 %
0066 sF = size(XF); 
0067 if sF(1) ~= nt | sF(2) ~= np,
0068     error('Input variables dimensions mismatch');
0069 end
0070 %
0071 % ----------------------------------------------------------------
0072 %
0073 XBT = abs(XBz);%Toroidal B-field
0074 XBP = sqrt(XBx.^2 + XBy.^2);%Poloidal B-field
0075 XB = sqrt(XBT.^2 + XBP.^2);%Total B-field
0076 %XB = sqrt(XBx.^2 + XBy.^2 + XBz.^2);
0077 %
0078 % calculation of the number of extrema of B
0079 %
0080 pnextr = sum(diff(sign(diff(XB([1:nt,2],:),1,1)),1,1) ~= 0,1);
0081 if pnextr ~= 2*ones(1,np),
0082     error('The magnetic field must have one minimum and one maximum only');
0083 end
0084 %
0085 % for each p point, the theta grid is temporarily shifted such that the minimum of B field
0086 % is located at i=1 (and also i=nt);
0087 %
0088 pBmin = (XB(1:nt-1,:) == repmat(min(XB,[],1),[nt-1,1]));
0089 pBmax = (XB(1:nt-1,:) == repmat(max(XB,[],1),[nt-1,1]));
0090 %
0091 if sum(pBmin,1) ~= ones(1,np),
0092     error('The magnetic field must have one minimum only');
0093 end
0094 if sum(pBmax,1) ~= ones(1,np),
0095     error('The magnetic field must have one maximum only');
0096 end
0097 %
0098 XCS1 = repmat((1:nt-1).',[1,np]);%cumsum(ones(nt-1,np),1);
0099 XCS2 = repmat(1:np,[nt-1,1]);%cumsum(ones(nt-1,np),2);
0100 %
0101 pit_Bmin = sum(pBmin.*XCS1);
0102 pit_Bmax = sum(pBmax.*XCS1);
0103 %
0104 % construction of the matrix of indices
0105 %
0106 Xit_Bmin = XCS1 + repmat(pit_Bmin,[nt-1,1]) - 1;
0107 Xit_Bmin(Xit_Bmin > nt-1) = Xit_Bmin(Xit_Bmin > nt-1) - (nt-1);
0108 %
0109 Xit_Bmin = Xit_Bmin + nt*(XCS2-1);
0110 %
0111 % permutation of the magnetic field matrices
0112 %
0113 XBx = XBx(Xit_Bmin);XBx = [XBx;XBx(1,:)];
0114 XBy = XBy(Xit_Bmin);XBy = [XBy;XBy(1,:)];
0115 XBz = XBz(Xit_Bmin);XBz = [XBz;XBz(1,:)];
0116 XB = XB(Xit_Bmin);XB = [XB;XB(1,:)];
0117 %
0118 XbBb = repmat(reshape(YBb,[1,np,nb]),[nt,1,1]);
0119 XbB = repmat(XB,[1,1,nb]);
0120 %
0121 XbBb = reshape(XbBb,[nt,np*nb]);
0122 XbB = reshape(XbB,[nt,np*nb]);
0123 %
0124 Ymask = XbB(1,:) <= XbBb(1,:) & sum(XbB > XbBb,1) > 0;
0125 Ymask1 = max(XbB,[],1) == XbBb(1,:);% particular case where B(theta) = max(XB)
0126 %
0127 XbCS = repmat((1:nt).',[1,np*nb]);%cumsum(ones(nt,np*nb));
0128 %
0129 Xbimax = (XbB > XbBb).*XbCS;
0130 YBb = XbBb(1,Ymask);
0131 nY = length(YBb);
0132 %
0133 % clear XbBb XbB XbCS
0134 %
0135 Xbimax(Xbimax == 0) = NaN;
0136 %
0137 Yi1max = min(Xbimax,[],1);
0138 Yi1min = Yi1max - 1;
0139 Yi2min = max(Xbimax,[],1);
0140 Yi2max = Yi2min + 1;
0141 %
0142 % clear Xbimax
0143 %
0144 %
0145 %xytheta1min =  theta(xyi1min);
0146 %xytheta1max =  theta(xyi1max);
0147 %xytheta2min =  theta(xyi2min);
0148 %xytheta2max =  theta(xyi2max);
0149 %
0150 Yipshift = nt*(repmat((1:np),[1,nb])-1);
0151 %
0152 YiX1min = Yi1min(Ymask) + Yipshift(Ymask);
0153 YiX1max = Yi1max(Ymask) + Yipshift(Ymask);
0154 YiX2min = Yi2min(Ymask) + Yipshift(Ymask);
0155 YiX2max = Yi2max(Ymask) + Yipshift(Ymask);
0156 %
0157 YBx1min = reshape(XBx(YiX1min),[1,nY]);
0158 YBx1max = reshape(XBx(YiX1max),[1,nY]);
0159 YBx2min = reshape(XBx(YiX2min),[1,nY]);
0160 YBx2max = reshape(XBx(YiX2max),[1,nY]);
0161 %
0162 YBy1min = reshape(XBy(YiX1min),[1,nY]);
0163 YBy1max = reshape(XBy(YiX1max),[1,nY]);
0164 YBy2min = reshape(XBy(YiX2min),[1,nY]);
0165 YBy2max = reshape(XBy(YiX2max),[1,nY]);
0166 %
0167 YBz1min = reshape(XBz(YiX1min),[1,nY]);
0168 YBz1max = reshape(XBz(YiX1max),[1,nY]);
0169 YBz2min = reshape(XBz(YiX2min),[1,nY]);
0170 YBz2max = reshape(XBz(YiX2max),[1,nY]);
0171 %
0172 YB1min = reshape(XB(YiX1min),[1,nY]);
0173 YB1max = reshape(XB(YiX1max),[1,nY]);
0174 YB2min = reshape(XB(YiX2min),[1,nY]);
0175 YB2max = reshape(XB(YiX2max),[1,nY]);
0176 %
0177 Y1 = YBx1min.*YBx1max + YBy1min.*YBy1max + YBz1min.*YBz1max;
0178 Y2 = YBx2min.*YBx2max + YBy2min.*YBy2max + YBz2min.*YBz2max;
0179 %
0180 Ydelta1 = sqrt((Y1 - YBb.^2).^2 - (YBb.^2 - YB1min.^2).*(YBb.^2 - YB1max.^2));
0181 Ydelta2 = sqrt((Y2 - YBb.^2).^2 - (YBb.^2 - YB2min.^2).*(YBb.^2 - YB2max.^2));
0182 %
0183 Yalpha1a = ((YB1min.^2 - Y1) + Ydelta1)./(YB1min.^2 + YB1max.^2 - 2*Y1);
0184 Yalpha1b = ((YB1min.^2 - Y1) - Ydelta1)./(YB1min.^2 + YB1max.^2 - 2*Y1);
0185 %
0186 Yalpha2a = ((YB2min.^2 - Y2) + Ydelta2)./(YB2min.^2 + YB2max.^2 - 2*Y2);
0187 Yalpha2b = ((YB2min.^2 - Y2) - Ydelta2)./(YB2min.^2 + YB2max.^2 - 2*Y2);
0188 %
0189 tol = 1e-5;
0190 %
0191 Ymask1a = Yalpha1a >= -tol & Yalpha1a <= 1+tol;
0192 Ymask1b = Yalpha1b >= -tol & Yalpha1b <= 1+tol;
0193 Ymask2a = Yalpha2a >= -tol & Yalpha2a <= 1+tol;
0194 Ymask2b = Yalpha2b >= -tol & Yalpha2b <= 1+tol;
0195 %
0196 if Ymask1a + Ymask1b ~= ones(1,length(Ymask1a)),
0197     error('root 1 not determined');
0198 end
0199 if Ymask2a + Ymask2b ~= ones(1,length(Ymask1a)),
0200     error('root 2 not determined');
0201 end
0202 %
0203 Yalpha1 = (Ymask1a.*Yalpha1a + ~Ymask1a.*Yalpha1b);
0204 Yalpha2 = (Ymask2a.*Yalpha2a + ~Ymask2a.*Yalpha2b);
0205 %
0206 Yit_Bmin = repmat(pit_Bmin,[1,nb]);
0207 Yit_Bmax = repmat(pit_Bmax,[1,nb]);
0208 %
0209 Ytheta1 = NaN*ones(np,nb);
0210 Ytheta2 = NaN*ones(np,nb);
0211 Ytheta1(Ymask1) = theta(Yit_Bmax(Ymask1));
0212 Ytheta2(Ymask1) = theta(Yit_Bmax(Ymask1));
0213 %
0214 Yi1min = Yi1min(Ymask) + Yit_Bmin(Ymask) - 1;
0215 Yi1max = Yi1max(Ymask) + Yit_Bmin(Ymask) - 1;
0216 Yi2min = Yi2min(Ymask) + Yit_Bmin(Ymask) - 1;
0217 Yi2max = Yi2max(Ymask) + Yit_Bmin(Ymask) - 1;
0218 %
0219 Yi1min(Yi1min > nt) = Yi1min(Yi1min > nt) - nt; 
0220 Yi1max(Yi1max > nt) = Yi1max(Yi1max > nt) - nt; 
0221 Yi2min(Yi2min > nt) = Yi2min(Yi2min > nt) - nt; 
0222 Yi2max(Yi2max > nt) = Yi2max(Yi2max > nt) - nt; 
0223 %
0224 Ytheta1(Ymask) = theta(Yi1min).'.*(1 - Yalpha1) + theta(Yi1max).'.*Yalpha1;
0225 Ytheta2(Ymask) = theta(Yi2min).'.*(1 - Yalpha2) + theta(Yi2max).'.*Yalpha2;
0226 %
0227 if nargout > 2,
0228     %
0229     if length(sF) == 2,
0230         sF = [sF,1];
0231     end
0232     %
0233     sFp = prod(sF(3:end));
0234     %
0235     YF1 = NaN*ones(np*nb,sFp);
0236     YF2 = NaN*ones(np*nb,sFp);
0237     %
0238     XF = reshape(XF,[nt*np,sFp]);
0239     %
0240     YF1(Ymask1,:) = XF(Yit_Bmax(Ymask1) + Yipshift(Ymask1),:);
0241     YF2(Ymask1,:) = XF(Yit_Bmax(Ymask1) + Yipshift(Ymask1),:);
0242     %
0243     Yalpha1 = repmat(Yalpha1.',[1,sFp]);
0244     Yalpha2 = repmat(Yalpha2.',[1,sFp]);    
0245     %
0246     YF1(Ymask,:) = XF(Yi1min + Yipshift(Ymask),:).*(1 - Yalpha1) + XF(Yi1max + Yipshift(Ymask),:).*Yalpha1;
0247     YF2(Ymask,:) = XF(Yi2min + Yipshift(Ymask),:).*(1 - Yalpha2) + XF(Yi2max + Yipshift(Ymask),:).*Yalpha2;
0248     %
0249     YF1 = reshape(YF1,[np,nb,sFp]);
0250     YF2 = reshape(YF2,[np,nb,sFp]);    
0251     %
0252 end
0253 
0254 
0255 
0256
interp_theta_jd

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
