testfitfluct_yp

PURPOSE ^

SYNOPSIS ^

function testfitfluct_yp(equil_fit,fluct,fluct_fit,rho,theta,mode)

DESCRIPTION ^

 Test the interpolated fluctuation spectrum

 INPUT
    - equil_fit: interpolated magnetic equilibrium structure
    - fluct: numerical plasma fluctuations structure
    - fluct_fit: interpolated plasma fluctuations structure
    - rho: normalized minor radius [m,1]
    - theta: poloidal angle (rad) [1,n]
    - mode: (1) -> (rho,theta), (2) -> (x,y) [1,1]
      (default = 1)

 OUTPUT: none

 By Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function testfitfluct_yp(equil_fit,fluct,fluct_fit,rho,theta,mode)
0002 %
0003 % Test the interpolated fluctuation spectrum
0004 %
0005 % INPUT
0006 %    - equil_fit: interpolated magnetic equilibrium structure
0007 %    - fluct: numerical plasma fluctuations structure
0008 %    - fluct_fit: interpolated plasma fluctuations structure
0009 %    - rho: normalized minor radius [m,1]
0010 %    - theta: poloidal angle (rad) [1,n]
0011 %    - mode: (1) -> (rho,theta), (2) -> (x,y) [1,1]
0012 %      (default = 1)
0013 %
0014 % OUTPUT: none
0015 %
0016 % By Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr)
0017 %
0018 if nargin < 4,
0019     rho = fluct.rho;
0020     dprho = gradient(fluct.rho(:));
0021     theta = fluct.theta;
0022     dtheta = gradient(fluct.theta);
0023 end
0024 %
0025 nrho = 1001;
0026 ntheta = 40001;
0027 %
0028 if nargin == 4, theta = [];end
0029 %
0030 if nargin < 6,mode = 1;end
0031 %
0032 if isempty(rho),rho = linspace(0,1,nrho);end
0033 if isempty(theta),theta = linspace(0,2*pi,ntheta);end
0034 %
0035 rho(abs(rho) > 1) = 1;%rho must be less than 1
0036 rho(rho < 0) = abs(rho(rho < 0));%rho must be positive
0037 %
0038 close all;
0039 %
0040 field_name = fieldnames(fluct_fit);
0041 %
0042 if (nargin >= 4) & (length(rho) == 1) & length(theta) > 1,%Calculation of the local magnetic equilibrium from equil_fit for the statistical analysis of the fluctuation spectrum
0043     x_a0_fit = ppval(equil_fit.x_fit.pp_a0,rho);
0044     x_an_fit = ppval(equil_fit.x_fit.pp_an,rho);
0045     x_bn_fit = ppval(equil_fit.x_fit.pp_bn,rho);
0046     x_da0drho_fit = ppval(equil_fit.x_fit.pp_da0drho,rho);
0047     x_dandrho_fit = ppval(equil_fit.x_fit.pp_dandrho,rho);
0048     x_dbndrho_fit = ppval(equil_fit.x_fit.pp_dbndrho,rho);
0049     %
0050     y_a0_fit = ppval(equil_fit.y_fit.pp_a0,rho);
0051     y_an_fit = ppval(equil_fit.y_fit.pp_an,rho);
0052     y_bn_fit = ppval(equil_fit.y_fit.pp_bn,rho);
0053     y_da0drho_fit = ppval(equil_fit.y_fit.pp_da0drho,rho);
0054     y_dandrho_fit = ppval(equil_fit.y_fit.pp_dandrho,rho);
0055     y_dbndrho_fit = ppval(equil_fit.y_fit.pp_dbndrho,rho);
0056     %
0057     Bx_a0_fit = ppval(equil_fit.Bx_fit.pp_a0,rho);
0058     Bx_an_fit = ppval(equil_fit.Bx_fit.pp_an,rho);
0059     Bx_bn_fit = ppval(equil_fit.Bx_fit.pp_bn,rho);
0060     Bx_da0drho_fit = ppval(equil_fit.Bx_fit.pp_da0drho,rho);
0061     Bx_dandrho_fit = ppval(equil_fit.Bx_fit.pp_dandrho,rho);
0062     Bx_dbndrho_fit = ppval(equil_fit.Bx_fit.pp_dbndrho,rho);
0063     %
0064     By_a0_fit = ppval(equil_fit.By_fit.pp_a0,rho);
0065     By_an_fit = ppval(equil_fit.By_fit.pp_an,rho);
0066     By_bn_fit = ppval(equil_fit.By_fit.pp_bn,rho);
0067     By_da0drho_fit = ppval(equil_fit.By_fit.pp_da0drho,rho);
0068     By_dandrho_fit = ppval(equil_fit.By_fit.pp_dandrho,rho);
0069     By_dbndrho_fit = ppval(equil_fit.By_fit.pp_dbndrho,rho);
0070     %
0071     BPHI_a0_fit = ppval(equil_fit.BPHI_fit.pp_a0,rho);
0072     BPHI_an_fit = ppval(equil_fit.BPHI_fit.pp_an,rho);
0073     BPHI_bn_fit = ppval(equil_fit.BPHI_fit.pp_bn,rho);
0074     BPHI_da0drho_fit = ppval(equil_fit.BPHI_fit.pp_da0drho,rho);
0075     BPHI_dandrho_fit = ppval(equil_fit.BPHI_fit.pp_dandrho,rho);
0076     BPHI_dbndrho_fit = ppval(equil_fit.BPHI_fit.pp_dbndrho,rho);
0077     %
0078     for ii = 1:length(theta), 
0079         [R_fit(:,ii),dRdtheta_fit(:,ii),d2Rdtheta2_fit(:,ii),dRdrho_fit(:,ii),d2Rdthetadrho_fit(:,ii)] = calcval_yp(equil_fit.x_fit,theta(ii),x_a0_fit(:),x_an_fit,x_bn_fit,x_da0drho_fit(:),x_dandrho_fit,x_dbndrho_fit);
0080         R_fit(:,ii) = R_fit(:,ii) + equil_fit.Rp;
0081         [Z_fit(:,ii),dZdtheta_fit(:,ii),d2Zdtheta2_fit(:,ii),dZdrho_fit(:,ii),d2Zdthetadrho_fit(:,ii)] = calcval_yp(equil_fit.y_fit,theta(ii),y_a0_fit(:),y_an_fit,y_bn_fit,y_da0drho_fit(:),y_dandrho_fit,y_dbndrho_fit);
0082         Z_fit(:,ii) = Z_fit(:,ii) + equil_fit.Zp;
0083         %
0084         r_fit(:,ii) = sqrt((R_fit(:,ii) - equil_fit.Rp).^2 + (Z_fit(:,ii) - equil_fit.Zp).^2);
0085         if (r_fit(1,ii) == 0), r_fit(1,ii) = prec; end
0086         drdtheta_fit(:,ii) = (dRdtheta_fit(:,ii).*(R_fit(:,ii) - equil_fit.Rp) + dZdtheta_fit(:,ii).*(Z_fit(:,ii) - equil_fit.Zp))./r_fit(:,ii);
0087         drdrho_fit(:,ii) = (dRdrho_fit(:,ii).*(R_fit(:,ii) - equil_fit.Rp) + dZdrho_fit(:,ii).*(Z_fit(:,ii) - equil_fit.Zp))./r_fit(:,ii);
0088         %
0089         salpha_fit(:,ii) = ((R_fit(:,ii) - equil_fit.Rp).*dRdtheta_fit(:,ii) + (Z_fit(:,ii) - equil_fit.Zp).*dZdtheta_fit(:,ii))./r_fit(:,ii)./sqrt(dRdtheta_fit(:,ii).^2 + dZdtheta_fit(:,ii).^2);
0090         dsalphadtheta_fit(:,ii) = ...
0091             sqrt(dRdtheta_fit(:,ii).^2 + dZdtheta_fit(:,ii).^2)./r_fit(:,ii) ...
0092             +((R_fit(:,ii) - equil_fit.Rp).*d2Rdtheta2_fit(:,ii) + (Z_fit(:,ii) - equil_fit.Zp).*d2Zdtheta2_fit(:,ii))./r_fit(:,ii)./sqrt(dRdtheta_fit(:,ii).^2 + dZdtheta_fit(:,ii).^2) ...
0093             -drdtheta_fit(:,ii).*((R_fit(:,ii) - equil_fit.Rp).*dRdtheta_fit(:,ii) + (Z_fit(:,ii) - equil_fit.Zp).*dZdtheta_fit(:,ii))./r_fit(:,ii).^2./sqrt(dRdtheta_fit(:,ii).^2 + dZdtheta_fit(:,ii).^2) ...
0094             -(dRdtheta_fit(:,ii).*d2Rdtheta2_fit(:,ii) + dZdtheta_fit(:,ii).*d2Zdtheta2_fit(:,ii)).*((R_fit(:,ii) - equil_fit.Rp).*dRdtheta_fit(:,ii) + (Z_fit(:,ii) - equil_fit.Zp).*dZdtheta_fit(:,ii))./r_fit(:,ii)./sqrt(dRdtheta_fit(:,ii).^2 + dZdtheta_fit(:,ii).^2)./(dRdtheta_fit(:,ii).^2 + dZdtheta_fit(:,ii).^2);     
0095         dsalphadrho_fit(:,ii) = ...
0096             (dRdtheta_fit(:,ii).*dRdrho_fit(:,ii) + dZdtheta_fit(:,ii).*dZdrho_fit(:,ii))./r_fit(:,ii)./sqrt(dRdtheta_fit(:,ii).^2 + dZdtheta_fit(:,ii).^2)...
0097             +((R_fit(:,ii) - equil_fit.Rp).*d2Rdthetadrho_fit(:,ii) + (Z_fit(:,ii) - equil_fit.Zp).*d2Zdthetadrho_fit(:,ii))./r_fit(:,ii)./sqrt(dRdtheta_fit(:,ii).^2 + dZdtheta_fit(:,ii).^2) ...
0098             -drdrho_fit(:,ii).*((R_fit(:,ii) - equil_fit.Rp).*dRdtheta_fit(:,ii) + (Z_fit(:,ii) - equil_fit.Zp).*dZdtheta_fit(:,ii))./r_fit(:,ii).^2./sqrt(dRdtheta_fit(:,ii).^2 + dZdtheta_fit(:,ii).^2) ...
0099             -(dRdtheta_fit(:,ii).*d2Rdthetadrho_fit(:,ii) + dZdtheta_fit(:,ii).*d2Zdthetadrho_fit(:,ii)).*((R_fit(:,ii) - equil_fit.Rp).*dRdtheta_fit(:,ii) + (Z_fit(:,ii) - equil_fit.Zp).*dZdtheta_fit(:,ii))./r_fit(:,ii)./sqrt(dRdtheta_fit(:,ii).^2 + dZdtheta_fit(:,ii).^2)./(dRdtheta_fit(:,ii).^2 + dZdtheta_fit(:,ii).^2);
0100         %
0101         calpha_fit(:,ii) = sqrt(1.0 - salpha_fit(:,ii).^2);
0102         if (calpha_fit(1,ii)== 0), calpha_fit(1,ii) = prec; end
0103         dcalphadtheta_fit(:,ii) = -salpha_fit(:,ii).*dsalphadtheta_fit(:,ii)./calpha_fit(:,ii);
0104         dcalphadrho_fit(:,ii) = -salpha_fit(:,ii).*dsalphadrho_fit(:,ii)./calpha_fit(:,ii);;
0105         %
0106         [BR_fit(:,ii),dBRdtheta_fit(:,ii),dummy,dBRdrho_fit(:,ii),dummy] = calcval_yp(equil_fit.Bx_fit,theta(ii),Bx_a0_fit(:),Bx_an_fit,Bx_bn_fit,Bx_da0drho_fit(:),Bx_dandrho_fit,Bx_dbndrho_fit);
0107         [BZ_fit(:,ii),dBZdtheta_fit(:,ii),dummy,dBZdrho_fit(:,ii),dummy] = calcval_yp(equil_fit.By_fit,theta(ii),By_a0_fit(:),By_an_fit,By_bn_fit,By_da0drho_fit(:),By_dandrho_fit,By_dbndrho_fit);
0108         [BPHI_fit(:,ii),dBPHIdtheta_fit(:,ii),dummy,dBPHIdrho_fit(:,ii),dummy] = calcval_yp(equil_fit.BPHI_fit,theta(ii),BPHI_a0_fit(:),BPHI_an_fit,BPHI_bn_fit,BPHI_da0drho_fit(:),BPHI_dandrho_fit,BPHI_dbndrho_fit);
0109         %
0110         BP_fit(:,ii) = sqrt(BR_fit(:,ii).^2 + BZ_fit(:,ii).^2);
0111         dBPdtheta_fit(:,ii) = (BR_fit(:,ii).*dBRdtheta_fit(:,ii) + BZ_fit(:,ii).*dBZdtheta_fit(:,ii))./BP_fit(:,ii);
0112         dBPdrho_fit(:,ii) = (BR_fit(:,ii).*dBRdrho_fit(:,ii) + BZ_fit(:,ii).*dBZdrho_fit(:,ii))./BP_fit(:,ii);
0113         %
0114         B_fit(:,ii) = sqrt(BR_fit(:,ii).^2 + BZ_fit(:,ii).^2 + BPHI_fit(:,ii).^2);
0115         dBdtheta_fit(:,ii) = (BR_fit(:,ii).*dBRdtheta_fit(:,ii) + BZ_fit(:,ii).*dBZdtheta_fit(:,ii) + BPHI_fit(:,ii).*dBPHIdtheta_fit(:,ii))./B_fit(:,ii);
0116         dBdrho_fit(:,ii) = (BR_fit(:,ii).*dBRdrho_fit(:,ii) + BZ_fit(:,ii).*dBZdrho_fit(:,ii) + BPHI_fit(:,ii).*dBPHIdrho_fit(:,ii))./B_fit(:,ii);
0117     end
0118     %
0119     P_fit = sign(equil_fit.psia)*BP_fit./B_fit;%With the sign of the plasma current (only valid for tokamak, not for RFP !)
0120     T_fit = BPHI_fit./B_fit;
0121 %
0122     a0_fit = ppval(fluct_fit.cm_fit.pp_a0,rho);
0123     an_fit = ppval(fluct_fit.cm_fit.pp_an,rho);
0124     bn_fit = ppval(fluct_fit.cm_fit.pp_bn,rho);
0125     da0drho_fit = ppval(fluct_fit.cm_fit.pp_da0drho,rho);
0126     dandrho_fit = ppval(fluct_fit.cm_fit.pp_dandrho,rho);
0127     dbndrho_fit = ppval(fluct_fit.cm_fit.pp_dbndrho,rho);
0128     %
0129     for ii = 1:length(theta),
0130         [cm_fit(:,ii)] = calcval_yp(fluct_fit.cm_fit,theta(ii),a0_fit(:),an_fit,bn_fit,da0drho_fit(:),dandrho_fit,dbndrho_fit);
0131         [int02theta_cm_fit(:,ii)] = intval_yp(fluct_fit.cm_fit,theta(ii),a0_fit(:),an_fit,bn_fit);%For the fluctuation phase calculation
0132     end
0133     %
0134     a0_fit = ppval(fluct_fit.cn_fit.pp_a0,rho);
0135     an_fit = ppval(fluct_fit.cn_fit.pp_an,rho);
0136     bn_fit = ppval(fluct_fit.cn_fit.pp_bn,rho);
0137     da0drho_fit = ppval(fluct_fit.cn_fit.pp_da0drho,rho);
0138     dandrho_fit = ppval(fluct_fit.cn_fit.pp_dandrho,rho);
0139     dbndrho_fit = ppval(fluct_fit.cn_fit.pp_dbndrho,rho);
0140     %
0141     for ii = 1:length(theta),
0142         [cn_fit(:,ii)] = calcval_yp(fluct_fit.cn_fit,theta(ii),a0_fit(:),an_fit,bn_fit,da0drho_fit(:),dandrho_fit,dbndrho_fit);
0143     end 
0144 %
0145     for i = 1:length(field_name),
0146         if isfield(fluct_fit,field_name(i)),
0147             if isstruct(fluct_fit.(char(field_name(i)))) & ~isempty(strfind(char(field_name(i)),'_fit')) & isempty(strfind(char(field_name(i)),'cm')) & isempty(strfind(char(field_name(i)),'cn')),
0148                 for ifield = 1:length(fluct_fit.(char(field_name(i)))),
0149                     subfield_name = fieldnames(fluct_fit.(char(field_name(i))));
0150                     v = regexprep(char(field_name(i)),'_fit','');
0151                     if fluct_fit.(char(field_name(i)))(ifield).model > 0,
0152                         kperp = fluct_fit.(char(field_name(i)))(ifield).kp;
0153                         nkperp = length(kperp);
0154                         phaseinit = fluct_fit.(char(field_name(i)))(ifield).phase;
0155                         %
0156                         j = 1;
0157                         for isubfield = 1:length(subfield_name),
0158                             if mode == 1,
0159                                 %
0160                                 % calculation at all radial positions of all coefficients for the Fourier series
0161                                 %
0162                                 if strfind(char(subfield_name(isubfield)),'_fit');
0163                                     a0_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_a0,rho);
0164                                     an_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_an,rho);
0165                                     bn_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_bn,rho);
0166                                     da0drho_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_da0drho,rho);
0167                                     dandrho_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_dandrho,rho);
0168                                     dbndrho_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_dbndrho,rho);
0169                                     %
0170                                     % calculation of the value at all poloidal angles
0171                                     %
0172                                     for ii = 1:length(theta),
0173                                         [w_fit{j}(:,ii)] = calcval_yp(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]),theta(ii),a0_fit(:),an_fit,bn_fit);
0174                                     end
0175                                     %
0176                                     j = j + 1;
0177                                 end
0178                             else
0179                                 error('Mode 2 (x,y) not yet implemented in fitfluct_yp.m');
0180                             end
0181                         end
0182                         %
0183                         if 20*max(kperp)*max(r_fit./T_fit./T_fit./sqrt(1 + calpha_fit.*calpha_fit.*P_fit.*P_fit./T_fit./T_fit))/equil_fit.ap > ntheta,
0184                             disp(['WARNING: The number of theta points is undestimated for a correct plot of the fluctuations. Please relaunch the test with ntheta > ',int2str(round(20*max(kperp)*max(r_fit./T_fit./T_fit./sqrt(1 + calpha_fit.*calpha_fit.*P_fit.*P_fit./T_fit./T_fit))/equil_fit.ap)),'.']);
0185                         end
0186                         %
0187                         dkperp = kperp(2) - kperp(1);%Uniform kperp grid assumed
0188                         %
0189                         if fluct_fit.(char(field_name(i)))(ifield).model == 1,%Gaussian model
0190                             %
0191                             phase = zeros(nkperp,length(theta));
0192                             %
0193                             phi = -theta.*P_fit.*r_fit./T_fit./R_fit./calpha_fit;%Phi is here chosen such that direction is perpendicular to the magnetic field line
0194                             %
0195                             for ikperp = 1:nkperp,
0196                                 phase(ikperp,:) = (int02theta_cm_fit + cn_fit(:,1).*phi)*kperp(ikperp) + phaseinit(ikperp);%toroidal angle phi is here chosen such that direction is perpendicular to the magnetic field line
0197                             end
0198                             %
0199                             if length(w_fit) == 3,
0200                                 df = zeros(nkperp,length(theta));
0201                                 %
0202                                 for ikperp = 1:nkperp,
0203                                     df(ikperp,:) = w_fit{1}.*cos(phase(ikperp,:)).*sqrt(2.*dkperp./w_fit{2}).*exp(-(kperp(ikperp) - w_fit{3}).^2./w_fit{2}./w_fit{2}./4)/(2*pi)^(0.25);
0204                                 end
0205                             else
0206                                 dfx = zeros(nkperp,length(theta));
0207                                 dfy = zeros(nkperp,length(theta));
0208                                 dfz = zeros(nkperp,length(theta));
0209                                 %
0210                                 for ikperp = 1:nkperp,
0211                                     dfx(ikperp,:) = w_fit{1}.*cos(phase(ikperp,:)).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25);
0212                                     dfy(ikperp,:) = w_fit{2}.*cos(phase(ikperp,:)).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25);
0213                                     dfz(ikperp,:) = w_fit{3}.*cos(phase(ikperp,:)).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25);
0214                                 end
0215                                 %
0216                                 df = sqrt(dfx.^2 + dfy.^2 + dfz.^2).*sign(dfz);
0217                             end
0218                             %
0219                             f = sum(df,1);
0220                             %
0221                             figure('Name',[v,'[',int2str(ifield),']']);
0222                             plot(theta,f,'r');
0223                             xlabel('\theta');ylabel([v,'[',int2str(ifield),']']);
0224                             title(['@ \rho = ',num2str(rho)]);
0225                             %
0226                             [h0,hx] = hist(f,max([11,round(ntheta/100)]));
0227                             h1 = h0/ntheta;
0228                             if length(w_fit) == 3,
0229                                 hm = (hx(2)-hx(1))*exp(-hx.*hx/2/w_fit{1}(1,1)^2)/sqrt(2*pi)/w_fit{1}(1,1);
0230                             else
0231                                 hm = (hx(2)-hx(1))*exp(-hx.*hx/2/(w_fit{1}(1,1)^2+w_fit{2}(1,1)^2+w_fit{3}(1,1)^2))/sqrt(2*pi)/sqrt((w_fit{1}(1,1)^2+w_fit{2}(1,1)^2+w_fit{3}(1,1)^2));
0232                             end
0233                             figure('Name',[v,'[',int2str(ifield),']']);
0234                             plot(hx,h1,'b',hx,hm,'r');
0235                             xlabel([v,'[',int2str(ifield),']']);
0236                             ylabel(['f(',v,'[',int2str(ifield),'])']);
0237                             title(['@ \rho = ',num2str(rho)]);
0238                             legend('Numerical distribution','1-D Maxwellian distribution');
0239                             drawnow;
0240                         end
0241                     else
0242                         disp('WARNING: no results for magnetic ripple.');
0243                     end    
0244                 end
0245             end
0246         end
0247     end
0248 else
0249     for i = 1:length(field_name)-1,
0250         if isfield(fluct_fit,field_name(i)),
0251             if isstruct(fluct_fit.(char(field_name(i)))) & ~isempty(strfind(char(field_name(i)),'_fit')),
0252                 for ifield = 1:length(fluct_fit.(char(field_name(i)))),
0253                     v = regexprep(char(field_name(i)),'_fit','');
0254                     subfield_name = fieldnames(fluct_fit.(char(field_name(i))));
0255                     if strfind(subfield_name{2},'pp_'),
0256                         if mode == 1,
0257                             %
0258                             % calculation at all radial positions of all coefficients for the Fourier series
0259                             %
0260                             a0_fit = ppval(fluct_fit.(char(field_name(i))).pp_a0,rho);
0261                             an_fit = ppval(fluct_fit.(char(field_name(i))).pp_an,rho);
0262                             bn_fit = ppval(fluct_fit.(char(field_name(i))).pp_bn,rho);
0263                             da0drho_fit = ppval(fluct_fit.(char(field_name(i))).pp_da0drho,rho);
0264                             dandrho_fit = ppval(fluct_fit.(char(field_name(i))).pp_dandrho,rho);
0265                             dbndrho_fit = ppval(fluct_fit.(char(field_name(i))).pp_dbndrho,rho);
0266                             %
0267                             % calculation of the value at all poloidal angles
0268                             %
0269                             for ii = 1:length(theta),
0270                                 [v_fit(:,ii),dvdtheta_fit(:,ii),d2vdtheta2_fit(:,ii),dvdrho_fit(:,ii),d2vdthetadrho_fit(:,ii)] = calcval_yp(fluct_fit.(char(field_name(i))),theta(ii),a0_fit(:),an_fit,bn_fit,da0drho_fit(:),dandrho_fit,dbndrho_fit);
0271                                 %
0272                                 if strfind(char(field_name(i)),'cm') & (length(rho) == 1) & length(theta) == 1,%for the fluctuation phase
0273                                     [iv_fit(:,ii),idvdrho_fit(:,ii)] = intval_yp(fluct_fit.(char(field_name(i))),theta(ii),a0_fit(:),an_fit,bn_fit,da0drho_fit(:),dandrho_fit,dbndrho_fit);
0274                                 end
0275                             end
0276                             %
0277                             if nargin < 4,
0278                                 %
0279                                 v0 = fluct.(regexprep(char(field_name(i)),'_fit',''));
0280                                 %
0281                                 [dvdtheta,dvdrho] = gradient(v0);
0282                                 dvdtheta = dvdtheta./(ones(length(dprho),1)*dtheta);
0283                                 dvdrho = dvdrho./(dprho(:)*ones(1,length(dtheta)));
0284                                 %
0285                                 [d2vdtheta2,d2vdthetadrho] = gradient(dvdtheta);
0286                                 d2vdtheta2 = d2vdtheta2./(ones(length(dprho),1)*dtheta);
0287                                 d2vdthetadrho = d2vdthetadrho./(dprho(:)*ones(1,length(dtheta)));
0288                                 %
0289                                 [d2vdrhodtheta,d2vdrho2] = gradient(dvdrho);
0290                                 d2vdrhodtheta = d2vdrhodtheta./(ones(length(dprho),1)*dtheta);
0291                                 d2vdrho2 = d2vdrho2./(dprho(:)*ones(1,length(dtheta)));
0292                                 %
0293                                 figure('Name',v);
0294                                 subplot(1,2,1),plot(rho,v0,'b+-',rho,v_fit,'r-');ylabel(v);xlabel('\rho');
0295                                 subplot(1,2,2),plot(theta/2/pi,v0,'b+-',theta/2/pi,v_fit,'r-');xlabel('\theta/2\pi');
0296                                 %
0297                                    figure('Name',['d',v,'/drho']);
0298                                 subplot(1,2,1),plot(rho,dvdrho,'b+-',rho,dvdrho_fit,'r-');ylabel(['d',v,'/drho']);xlabel('\rho');
0299                                 subplot(1,2,2),plot(theta/2/pi,dvdrho,'b+-',theta/2/pi,dvdrho_fit,'r-');xlabel('\theta/2\pi');
0300                                 %
0301                                 figure('Name',['d',v,'/dtheta']);
0302                                 subplot(1,2,1),plot(rho,dvdtheta,'b+-',rho,dvdtheta_fit,'r-');ylabel(['d',v,'/dtheta']);xlabel('\rho');
0303                                 subplot(1,2,2),plot(theta/2/pi,dvdtheta,'b+-',theta/2/pi,dvdtheta_fit,'r-');xlabel('\theta/2\pi');
0304                                 %
0305                                 figure('Name',['d2',v,'/dtheta2']);
0306                                 subplot(1,2,1),plot(rho,d2vdtheta2,'b+-',rho,d2vdtheta2_fit,'r-');ylabel(['d2',v,'/dtheta2']);xlabel('\rho');
0307                                 subplot(1,2,2),plot(theta/2/pi,d2vdtheta2,'b+-',theta/2/pi,d2vdtheta2_fit,'r-');xlabel('\theta/2\pi');
0308                                 %
0309                                 figure('Name',['d2',v,'/dthetadrho']);
0310                                 subplot(1,2,1),plot(rho,d2vdthetadrho,'b+-',rho,d2vdthetadrho_fit,'r-');ylabel(['d2',v,'/dthetadrho']);xlabel('\rho');
0311                                 subplot(1,2,2),plot(theta/2/pi,d2vdthetadrho,'b+-',theta/2/pi,d2vdthetadrho_fit,'r-');xlabel('\theta/2\pi');
0312                             else
0313                                 if (length(rho) > 1) & (length(theta) > 1),
0314                                     figure('Name',[v,'[',int2str(ifield),'] - ',char(field_name(i))]);
0315                                     subplot(1,2,1),plot(rho,v_fit,'r-');ylabel([char(field_name(i))]);xlabel('\rho');
0316                                     subplot(1,2,2),plot(theta/2/pi,v_fit,'r-');xlabel('\theta/2\pi');
0317                                     %
0318                                     figure('Name',[v,'[',int2str(ifield),']: d',char(field_name(i)),'/drho']);
0319                                     subplot(1,2,1),plot(rho,dvdrho_fit,'r-');ylabel(['d',char(subfield_name(isubfield)),'/drho']);xlabel('\rho');
0320                                     subplot(1,2,2),plot(theta/2/pi,dvdrho_fit,'r-');xlabel('\theta/2\pi');
0321                                     %
0322                                     figure('Name',[v,'[',int2str(ifield),'] d',char(field_name(i)),'/dtheta']);
0323                                     subplot(1,2,1),plot(rho,dvdtheta_fit,'r-');ylabel(['d',char(subfield_name(isubfield)),'/dtheta']);xlabel('\rho');
0324                                     subplot(1,2,2),plot(theta/2/pi,dvdtheta_fit,'r-');xlabel('\theta/2\pi');
0325                                     %
0326                                     figure('Name',[v,'[',int2str(ifield),'] d2',char(field_name(i)),'/dtheta2']);
0327                                     subplot(1,2,1),plot(rho,d2vdtheta2_fit,'r-');ylabel(['d2',char(subfield_name(isubfield)),'/dtheta2']);xlabel('\rho');
0328                                     subplot(1,2,2),plot(theta/2/pi,d2vdtheta2_fit,'r-');xlabel('\theta/2\pi');
0329                                     %
0330                                     figure('Name',[v,'[',int2str(ifield),'] d2',char(field_name(i)),'/dthetadrho']);
0331                                     subplot(1,2,1),plot(rho,d2vdthetadrho_fit,'r-');ylabel(['d2',char(subfield_name(isubfield)),'/dthetadrho']);xlabel('\rho');
0332                                     subplot(1,2,2),plot(theta/2/pi,d2vdthetadrho_fit,'r-');xlabel('\theta/2\pi');
0333                                     %
0334                                 elseif (length(rho) == 1) & length(theta) == 1,
0335                                     disp(['-------------- @ (rho = ',num2str(rho),',theta = ',num2str(theta),') --------------']);
0336     
0337                                     disp([char(field_name(i)),' = ',num2str(v_fit)]);
0338                                     disp(['d',char(field_name(i)),'/drho = ',num2str(dvdrho_fit)]);
0339                                     disp(['d',char(field_name(i)),'/dtheta = ',num2str(dvdtheta_fit)]);
0340                                     % disp(['d2',char(field_name(i)),'/drho2 = ',num2str()]);
0341                                     disp(['d2',char(field_name(i)),'/dtheta2 = ',num2str(d2vdtheta2_fit)]);
0342                                     disp(['d2',char(field_name(i)),'/dthetadrho = ',num2str(d2vdthetadrho_fit)]);
0343                                     %
0344                                     if strfind(char(field_name(i)),'cm'),%for the fluctuation phase
0345                                         disp(['i',char(field_name(i)),' = ',num2str(iv_fit)]);
0346                                         disp(['id',char(field_name(i)),'/drho = ',num2str(idvdrho_fit)]);
0347                                     end
0348                                 end
0349                             end    
0350                         else
0351                             error('Mode 2 (x,y) not yet implemented in fitfluct_yp.m');
0352                         end
0353                     else
0354                         for isubfield = 3:length(subfield_name)-1,
0355                             if mode == 1,
0356                                 %
0357                                 % calculation at all radial positions of all coefficients for the Fourier series
0358                                 %
0359                                 a0_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_a0,rho);
0360                                 an_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_an,rho);
0361                                 bn_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_bn,rho);
0362                                 da0drho_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_da0drho,rho);
0363                                 dandrho_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_dandrho,rho);
0364                                 dbndrho_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_dbndrho,rho);
0365                                 %
0366                                 % calculation of the value at all poloidal angles
0367                                 %
0368                                 for ii = 1:length(theta),
0369                                     [v_fit(:,ii),dvdtheta_fit(:,ii),d2vdtheta2_fit(:,ii),dvdrho_fit(:,ii),d2vdthetadrho_fit(:,ii)] = calcval_yp(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]),theta(ii),a0_fit(:),an_fit,bn_fit,da0drho_fit(:),dandrho_fit,dbndrho_fit);
0370                                 end
0371                                 %
0372                                 if nargin < 4,
0373                                     %
0374                                     v0 = fluct.(regexprep(char(field_name(i)),'_fit',''))(ifield).(regexprep(char(subfield_name(isubfield)),'_fit',''));
0375                                     %
0376                                     [dvdtheta,dvdrho] = gradient(v0);
0377                                     dvdtheta = dvdtheta./(ones(length(dprho),1)*dtheta);
0378                                     dvdrho = dvdrho./(dprho(:)*ones(1,length(dtheta)));
0379                                     %
0380                                     [d2vdtheta2,d2vdthetadrho] = gradient(dvdtheta);
0381                                     d2vdtheta2 = d2vdtheta2./(ones(length(dprho),1)*dtheta);
0382                                     d2vdthetadrho = d2vdthetadrho./(dprho(:)*ones(1,length(dtheta)));
0383                                     %
0384                                     [d2vdrhodtheta,d2vdrho2] = gradient(dvdrho);
0385                                     d2vdrhodtheta = d2vdrhodtheta./(ones(length(dprho),1)*dtheta);
0386                                     d2vdrho2 = d2vdrho2./(dprho(:)*ones(1,length(dtheta)));
0387                                     %
0388                                     figure('Name',[v,'[',int2str(ifield),'] - ',char(subfield_name(isubfield))]);
0389                                     subplot(1,2,1),plot(rho,v0,'b+-',rho,v_fit,'r-');ylabel([char(subfield_name(isubfield))]);xlabel('\rho');
0390                                     subplot(1,2,2),plot(theta/2/pi,v0,'b+-',theta/2/pi,v_fit,'r-');xlabel('\theta/2\pi');
0391                                     %
0392                                        figure('Name',[v,'[',int2str(ifield),']: d',char(subfield_name(isubfield)),'/drho']);
0393                                     subplot(1,2,1),plot(rho,dvdrho,'b+-',rho,dvdrho_fit,'r-');ylabel(['d',char(subfield_name(isubfield)),'/drho']);xlabel('\rho');
0394                                     subplot(1,2,2),plot(theta/2/pi,dvdrho,'b+-',theta/2/pi,dvdrho_fit,'r-');xlabel('\theta/2\pi');
0395                                     %
0396                                     figure('Name',[v,'[',int2str(ifield),'] d',char(subfield_name(isubfield)),'/dtheta']);
0397                                     subplot(1,2,1),plot(rho,dvdtheta,'b+-',rho,dvdtheta_fit,'r-');ylabel(['d',char(subfield_name(isubfield)),'/dtheta']);xlabel('\rho');
0398                                     subplot(1,2,2),plot(theta/2/pi,dvdtheta,'b+-',theta/2/pi,dvdtheta_fit,'r-');xlabel('\theta/2\pi');
0399                                     %
0400                                     figure('Name',[v,'[',int2str(ifield),'] d2',char(subfield_name(isubfield)),'/dtheta2']);
0401                                     subplot(1,2,1),plot(rho,d2vdtheta2,'b+-',rho,d2vdtheta2_fit,'r-');ylabel(['d2',char(subfield_name(isubfield)),'/dtheta2']);xlabel('\rho');
0402                                     subplot(1,2,2),plot(theta/2/pi,d2vdtheta2,'b+-',theta/2/pi,d2vdtheta2_fit,'r-');xlabel('\theta/2\pi');
0403                                     %
0404                                     figure('Name',[v,'[',int2str(ifield),'] d2',char(subfield_name(isubfield)),'/dthetadrho']);
0405                                     subplot(1,2,1),plot(rho,d2vdthetadrho,'b+-',rho,d2vdthetadrho_fit,'r-');ylabel(['d2',char(subfield_name(isubfield)),'/dthetadrho']);xlabel('\rho');
0406                                     subplot(1,2,2),plot(theta/2/pi,d2vdthetadrho,'b+-',theta/2/pi,d2vdthetadrho_fit,'r-');xlabel('\theta/2\pi');
0407                                 else
0408                                     if (length(rho) > 1) & (length(theta) > 1),
0409                                         figure('Name',[v,'[',int2str(ifield),'] - ',char(subfield_name(isubfield))]);
0410                                         subplot(1,2,1),plot(rho,v_fit,'r-');ylabel([char(subfield_name(isubfield))]);xlabel('\rho');
0411                                         subplot(1,2,2),plot(theta/2/pi,v_fit,'r-');xlabel('\theta/2\pi');
0412                                         %
0413                                         figure('Name',[v,'[',int2str(ifield),']: d',char(subfield_name(isubfield)),'/drho']);
0414                                         subplot(1,2,1),plot(rho,dvdrho_fit,'r-');ylabel(['d',char(subfield_name(isubfield)),'/drho']);xlabel('\rho');
0415                                         subplot(1,2,2),plot(theta/2/pi,dvdrho_fit,'r-');xlabel('\theta/2\pi');
0416                                         %
0417                                         figure('Name',[v,'[',int2str(ifield),'] d',char(subfield_name(isubfield)),'/dtheta']);
0418                                         subplot(1,2,1),plot(rho,dvdtheta_fit,'r-');ylabel(['d',char(subfield_name(isubfield)),'/dtheta']);xlabel('\rho');
0419                                         subplot(1,2,2),plot(theta/2/pi,dvdtheta_fit,'r-');xlabel('\theta/2\pi');
0420                                         %
0421                                         figure('Name',[v,'[',int2str(ifield),'] d2',char(subfield_name(isubfield)),'/dtheta2']);
0422                                         subplot(1,2,1),plot(rho,d2vdtheta2_fit,'r-');ylabel(['d2',char(subfield_name(isubfield)),'/dtheta2']);xlabel('\rho');
0423                                         subplot(1,2,2),plot(theta/2/pi,d2vdtheta2_fit,'r-');xlabel('\theta/2\pi');
0424                                         %
0425                                         figure('Name',[v,'[',int2str(ifield),'] d2',char(subfield_name(isubfield)),'/dthetadrho']);
0426                                         subplot(1,2,1),plot(rho,d2vdthetadrho_fit,'r-');ylabel(['d2',char(subfield_name(isubfield)),'/dthetadrho']);xlabel('\rho');
0427                                         subplot(1,2,2),plot(theta/2/pi,d2vdthetadrho_fit,'r-');xlabel('\theta/2\pi');
0428                                         %
0429                                     elseif (length(rho) == 1) & length(theta) == 1,
0430                                         disp(['-------------- ',v,'[',int2str(ifield),']  @ (rho = ',num2str(rho),',theta = ',num2str(theta),') --------------']);
0431                                         disp([char(subfield_name(isubfield)),'  = ',num2str(v_fit)]);
0432                                         disp(['d',char(subfield_name(isubfield)),'/drho = ',num2str(dvdrho_fit)]);
0433                                         disp(['d',char(subfield_name(isubfield)),'/dtheta = ',num2str(dvdtheta_fit)]);
0434                                         % disp(['d2',char(subfield_name(isubfield)),'/drho2 = ',num2str()]);
0435                                         disp(['d2',char(subfield_name(isubfield)),'/dtheta2 = ',num2str(d2vdtheta2_fit)]);
0436                                         disp(['d2',char(subfield_name(isubfield)),'/dthetadrho = ',num2str(d2vdthetadrho_fit)]);
0437                                     end
0438                                 end
0439                             else
0440                                 error('Mode 2 (x,y) not yet implemented in fitfluct_yp.m');
0441                             end
0442                         end
0443                     end
0444                 end
0445             end                
0446         end                
0447     end
0448 end
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.