fluctval_yp

PURPOSE ^

SYNOPSIS ^

function [fluctval] = fluctval_yp(equil_fit,fluct_fit,rho,theta,phi,theta0,phi0)

DESCRIPTION ^

 This function calculates the local perturbation of the magnetic equilibrium at any position in the plasma

   INPUT :

       - equil_fit: fitted axisymmetric magnetic equilibrium structure
       - fluct_fit: structure of the magnetic equilibrium fluctuations (including magnetic ripple)
       - rho: normalized minor radius (a.u.) [1,1]
       - theta: poloidal angle (rad) [1,1]
       - phi: toroidal angle (rad) [1,1]
       - theta0: reference of the initial poloidal angle for some fluctuation models (rad) [1,1]
       - phi0: reference of the initial toroidal angle for some fluctuation models (rad) [1,1]

   OUTPUT :

       - fluctval: structure of the values of the magnetic equilibrium fluctuations at (rho, theta, phi)

 by Y. Peysson DRFC-CEA <yves.peysson@cea.fr> and J. Decker DRFC-CEA <joan.decker@cea.fr>

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [fluctval] = fluctval_yp(equil_fit,fluct_fit,rho,theta,phi,theta0,phi0)
0002 %
0003 % This function calculates the local perturbation of the magnetic equilibrium at any position in the plasma
0004 %
0005 %   INPUT :
0006 %
0007 %       - equil_fit: fitted axisymmetric magnetic equilibrium structure
0008 %       - fluct_fit: structure of the magnetic equilibrium fluctuations (including magnetic ripple)
0009 %       - rho: normalized minor radius (a.u.) [1,1]
0010 %       - theta: poloidal angle (rad) [1,1]
0011 %       - phi: toroidal angle (rad) [1,1]
0012 %       - theta0: reference of the initial poloidal angle for some fluctuation models (rad) [1,1]
0013 %       - phi0: reference of the initial toroidal angle for some fluctuation models (rad) [1,1]
0014 %
0015 %   OUTPUT :
0016 %
0017 %       - fluctval: structure of the values of the magnetic equilibrium fluctuations at (rho, theta, phi)
0018 %
0019 % by Y. Peysson DRFC-CEA <yves.peysson@cea.fr> and J. Decker DRFC-CEA <joan.decker@cea.fr>
0020 %
0021 if nargin < 5,
0022     error('Not enough input arguments in fluctval_yp.m')
0023 end
0024 %
0025 if nargin == 5,
0026     theta0 = theta;
0027     phi0 = phi;
0028 end
0029 %
0030 if nargin == 6,
0031     phi0 = phi;
0032 end
0033 %
0034 equilval = equilval_yp(equil_fit,rho,theta);%axisymmetric magnetic equilibrium
0035 %
0036 Te = equilval.Te;
0037 dTedrho = equilval.dTedrho;
0038 ne = equilval.ne;
0039 dnedrho = equilval.dnedrho;
0040 zTi = equilval.zTi;
0041 dzTidrho = equilval.dzTidrho;
0042 zni = equilval.zni;
0043 dznidrho = equilval.dznidrho;
0044 R = equilval.Rp + equilval.x;
0045 r = equilval.r;
0046 salpha = equilval.salpha;
0047 dsalphadtheta = equilval.dsalphadtheta;
0048 dsalphadrho = equilval.dsalphadrho;
0049 calpha = equilval.calpha;
0050 dcalphadtheta = equilval.dcalphadtheta;
0051 dcalphadrho = equilval.dcalphadrho;
0052 BPHI = equilval.Bz;
0053 dBPHIdtheta = equilval.dBzdtheta;
0054 dBPHIdrho = equilval.dBzdrho;
0055 BP = equilval.BP;
0056 B = equilval.B;
0057 dBdtheta = equilval.dBdtheta;
0058 dBdrho = equilval.dBdrho;
0059 gradrho = equilval.gradrho;
0060 dgradrhodtheta = equilval.dgradrhodtheta;
0061 dgradrhodrho = equilval.dgradrhodrho;
0062 %
0063 fluctval.ner = 0;
0064 fluctval.dnerdrho = 0;
0065 fluctval.dnerdtheta = 0;
0066 fluctval.dnerdphi = 0;
0067 %
0068 fluctval.Ter = 0;
0069 fluctval.dTerdrho = 0;
0070 fluctval.dTerdtheta = 0;
0071 fluctval.dTerdphi = 0;
0072 %
0073 fluctval.znir = zeros(size(zni'));
0074 fluctval.dznirdrho = zeros(size(zni'));
0075 fluctval.dznirdtheta = zeros(size(zni'));
0076 fluctval.dznirdphi = zeros(size(zni'));
0077 %
0078 fluctval.zTir = zeros(size(zTi'));
0079 fluctval.dzTirdrho = zeros(size(zTi'));
0080 fluctval.dzTirdtheta = zeros(size(zTi'));
0081 fluctval.dzTirdphi = zeros(size(zTi'));
0082 %
0083 fluctval.Bxr = 0;
0084 fluctval.dBxrdrho = 0;
0085 fluctval.dBxrdtheta = 0;
0086 fluctval.dBxrdphi = 0;
0087 %
0088 fluctval.Byr = 0;
0089 fluctval.dByrdrho = 0;
0090 fluctval.dByrdtheta = 0;
0091 fluctval.dByrdphi = 0;
0092 %
0093 fluctval.Bzr = 0;
0094 fluctval.dBzrdrho = 0;
0095 fluctval.dBzrdtheta = 0;
0096 fluctval.dBzrdphi = 0;
0097 %
0098 field_name = fieldnames(fluct_fit);
0099 %
0100 L_rho = fluct_fit.L_rho;
0101 dL_rhodrho = 0;
0102 dL_rhodtheta = 0;
0103 dL_rhodphi = 0;
0104 %
0105 xq_fit = ppval(fluct_fit.xq_fit.pp_f,rho);
0106 dxqdrho_fit = ppval(fluct_fit.xq_fit.pp_dfdrho,rho);
0107 dxqdtheta_fit = 0;
0108 dxqdphi_fit = 0;
0109 %
0110 xL_theta_fit = ppval(fluct_fit.xL_theta_fit.pp_f,rho);
0111 dxL_thetadrho_fit = ppval(fluct_fit.xL_theta_fit.pp_dfdrho,rho);
0112 dxL_thetadtheta_fit = 0;
0113 dxL_thetadphi_fit = 0; 
0114 %
0115 xL_perp_fit = ppval(fluct_fit.xL_perp_fit.pp_f,rho);
0116 dxL_perpdrho_fit = ppval(fluct_fit.xL_perp_fit.pp_dfdrho,rho);
0117 dxL_perpdtheta_fit = 0;
0118 dxL_perpdphi_fit = 0;
0119 %
0120 a0_fit = ppval(fluct_fit.XL_phi_fit.pp_a0,rho);
0121 an_fit = ppval(fluct_fit.XL_phi_fit.pp_an,rho);
0122 bn_fit = ppval(fluct_fit.XL_phi_fit.pp_bn,rho);
0123 da0drho_fit = ppval(fluct_fit.XL_phi_fit.pp_da0drho,rho);
0124 dandrho_fit = ppval(fluct_fit.XL_phi_fit.pp_dandrho,rho);
0125 dbndrho_fit = ppval(fluct_fit.XL_phi_fit.pp_dbndrho,rho);
0126 %
0127 for ii = 1:length(theta),
0128     [XL_phi_fit(:,ii),dXL_phidtheta_fit(:,ii),dummy,dXL_phidrho_fit(:,ii)] = calcval_yp(fluct_fit.XL_phi_fit,theta(ii),a0_fit(:),an_fit,bn_fit,da0drho_fit(:),dandrho_fit,dbndrho_fit);
0129 end
0130 dXL_phidphi_fit = 0;
0131 %
0132 a0_fit = ppval(fluct_fit.Xcm_fit.pp_a0,rho);
0133 an_fit = ppval(fluct_fit.Xcm_fit.pp_an,rho);
0134 bn_fit = ppval(fluct_fit.Xcm_fit.pp_bn,rho);
0135 da0drho_fit = ppval(fluct_fit.Xcm_fit.pp_da0drho,rho);
0136 dandrho_fit = ppval(fluct_fit.Xcm_fit.pp_dandrho,rho);
0137 dbndrho_fit = ppval(fluct_fit.Xcm_fit.pp_dbndrho,rho);
0138 %
0139 for ii = 1:length(theta),
0140     [Xcm_fit(:,ii),dXcmdtheta_fit(:,ii),dummy,dXcmdrho_fit(:,ii)] = calcval_yp(fluct_fit.Xcm_fit,theta(ii),a0_fit(:),an_fit,bn_fit,da0drho_fit(:),dandrho_fit,dbndrho_fit);%For the fluctuation phase calculation (model 5)
0141     [Xicm_fit(:,ii),dXicmdrho_fit(:,ii)] = intval_yp(fluct_fit.Xcm_fit,theta0,theta(ii),a0_fit(:),an_fit,bn_fit,da0drho_fit(:),dandrho_fit,dbndrho_fit);%For the fluctuation phase calculation (model 5)
0142 end
0143 dXcmdphi_fit = 0;
0144 %
0145 a0_fit = ppval(fluct_fit.Xcn_fit.pp_a0,rho);
0146 an_fit = ppval(fluct_fit.Xcn_fit.pp_an,rho);
0147 bn_fit = ppval(fluct_fit.Xcn_fit.pp_bn,rho);
0148 da0drho_fit = ppval(fluct_fit.Xcn_fit.pp_da0drho,rho);
0149 dandrho_fit = ppval(fluct_fit.Xcn_fit.pp_dandrho,rho);
0150 dbndrho_fit = ppval(fluct_fit.Xcn_fit.pp_dbndrho,rho);
0151 %
0152 [Xcn0_fit(:,ii),dXcn0dtheta_fit(:,ii),dummy,dXcn0drho_fit(:,ii)] = calcval_yp(fluct_fit.Xcn_fit,theta0,a0_fit(:),an_fit,bn_fit,da0drho_fit(:),dandrho_fit,dbndrho_fit);%For the fluctuation phase calculation (model 5)
0153 dXcn0dphi_fit = 0;
0154 %
0155 for i = 1:length(field_name),
0156     if isfield(fluct_fit,field_name(i)),
0157         if isstruct(fluct_fit.(char(field_name(i)))) && (strcmp(char(field_name(i)),'B_fit') || strcmp(char(field_name(i)),'ne_fit')),   
0158             for ifield = 1:length(fluct_fit.(char(field_name(i)))),
0159                 subfield_name = fieldnames(fluct_fit.(char(field_name(i))));
0160                 v = regexprep(char(field_name(i)),'_fit','');
0161                 %
0162                 sigmar_max = fluct_fit.(char(field_name(i)))(ifield).sigmar_max;
0163                 sigmar_hwhm = fluct_fit.(char(field_name(i)))(ifield).sigmar_hwhm;
0164                 sigmar_rho = fluct_fit.(char(field_name(i)))(ifield).sigmar_rho;
0165                 polmode = fluct_fit.(char(field_name(i)))(ifield).polmode;
0166                 epsi_rho = fluct_fit.(char(field_name(i)))(ifield).epsi_rho;
0167                 epsi_theta = fluct_fit.(char(field_name(i)))(ifield).epsi_theta;
0168                 epsi_phi = fluct_fit.(char(field_name(i)))(ifield).epsi_phi;
0169                 lmin = fluct_fit.(char(field_name(i)))(ifield).lmin;
0170                 mmin = fluct_fit.(char(field_name(i)))(ifield).mmin;
0171                 nmin = fluct_fit.(char(field_name(i)))(ifield).nmin;
0172                 lmax = fluct_fit.(char(field_name(i)))(ifield).lmax;
0173                 mmax = fluct_fit.(char(field_name(i)))(ifield).mmax;
0174                 nmax = fluct_fit.(char(field_name(i)))(ifield).nmax;            
0175                 phaseinit = fluct_fit.(char(field_name(i)))(ifield).phase;
0176                 %
0177                 j = 1;
0178                 for isubfield = 1:length(subfield_name),
0179                     %
0180                     % calculation at all radial positions of all coefficients for the Fourier series
0181                     %
0182                     if strfind(char(subfield_name(isubfield)),'_fit');
0183                         a0_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_a0,rho);
0184                         an_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_an,rho);
0185                         bn_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_bn,rho);
0186                         da0drho_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_da0drho,rho);
0187                         dandrho_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_dandrho,rho);
0188                         dbndrho_fit = ppval(fluct_fit.(char(field_name(i)))(ifield).([char(subfield_name(isubfield))]).pp_dbndrho,rho);
0189                         %
0190                         % calculation of the value at all poloidal angles
0191                         %
0192                         for ii = 1:length(theta),
0193                             [w_fit{j}(:,ii),dwdtheta_fit{j}(:,ii),d2wdtheta2_fit{j}(:,ii),dwdrho_fit{j}(:,ii),d2wdthetadrho_fit{j}(:,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);
0194                         end
0195                         %
0196                         j = j + 1;
0197                     end
0198                 end
0199                 %
0200                 if fluct_fit.(char(field_name(i)))(ifield).model == 7,%(1) -> 2-D Gaussian drift-wave like model local (rho,lperp), absolute epsi values (m)
0201                     if strcmp(char(field_name(i)),'ne_fit'),
0202 
0203                     elseif strcmp(char(field_name(i)),'B_fit'),
0204                         
0205                         
0206                         
0207                         
0208                     end
0209                 elseif fluct_fit.(char(field_name(i)))(ifield).model == 6,%(1) -> 2-D Gaussian drift-wave like model local (rho,cn,cm), absolute epsi values (m)
0210                     if strcmp(char(field_name(i)),'ne_fit'),
0211                         %
0212                         [l,m] = build_tensor_yp(lmin:lmax,mmin:mmax);
0213                         %
0214                         u_l = pi*epsi_rho*[lmin:lmax]/L_rho;
0215                         u_m = pi*epsi_theta*[mmin:mmax]/xL_perp_fit;%here, epsi_perp = epsi_theta
0216                         %
0217                         du_ldrho = -u_l/L_rho*dL_rhodrho;
0218                         du_mdrho = -u_m/xL_perp_fit*dxL_perpdrho_fit;
0219                         %
0220                         du_ldtheta = -u_l/L_rho*dL_rhodtheta;
0221                         du_mdtheta = -u_m/xL_perp_fit*dxL_perpdtheta_fit;
0222                         %
0223                         du_ldphi = -u_l/L_rho*dL_rhodphi;
0224                         du_mdphi = -u_m/xL_perp_fit*dxL_perpdphi_fit;
0225                         %
0226                         expu_l = exp(-u_l.^2/2);
0227                         expu_m = exp(-u_m.^2/2);
0228                         %
0229                         [al,am] = build_tensor_yp(expu_l,expu_m);
0230                         alm = al.*am;
0231                         %
0232                         dexpu_ldrho = -expu_l.*u_l.*du_ldrho;
0233                         dexpu_mdrho = -expu_m.*u_m.*du_mdrho;
0234                         %
0235                         [daldrho,damdrho] = build_tensor_yp(dexpu_ldrho,dexpu_mdrho);
0236                         dalmdrho = daldrho.*am + al.*damdrho;
0237                         %
0238                         dexpu_ldtheta = -expu_l.*u_l.*du_ldtheta;
0239                         dexpu_mdtheta = -expu_m.*u_m.*du_mdtheta;
0240                         %
0241                         [daldtheta,damdtheta] = build_tensor_yp(dexpu_ldtheta,dexpu_mdtheta);
0242                         dalmdtheta = daldtheta.*am + al.*damdtheta;                        
0243                         %
0244                         dexpu_ldphi = -expu_l.*u_l.*du_ldphi;
0245                         dexpu_mdphi = -expu_m.*u_m.*du_mdphi;
0246                         %
0247                         [daldphi,damdphi] = build_tensor_yp(dexpu_ldphi,dexpu_mdphi);
0248                         dalmdphi = daldphi.*am + al.*damdphi;
0249                         %
0250                         kperp = 2*pi*m./xL_perp_fit;
0251                         dkperpdrho = -kperp./xL_perp_fit.*dxL_perpdrho_fit;
0252                         dkperpdtheta = 0;
0253                         dkperpdphi = 0;
0254                         %
0255                         phase = 2*pi*l*rho + (Xicm_fit + Xcn0_fit*(phi-phi0))*kperp + reshape(phaseinit,length(lmin:lmax),length(mmin:mmax));
0256                         dphasedrho = 2*pi*l + (dXicmdrho_fit + dXcn0drho_fit*(phi-phi0))*kperp + (Xicm_fit + Xcn0_fit*(phi-phi0))*dkperpdrho;
0257                         dphasedtheta = Xcm_fit*kperp + (Xicm_fit + Xcn0_fit*(phi-phi0))*dkperpdtheta;
0258                         dphasedphi = Xcn0_fit*kperp + (Xicm_fit + Xcn0_fit*(phi-phi0))*dkperpdphi;
0259                         %
0260                         Flambda = (1.0 + polmode + (1.0 - polmode)*cos(theta))/2.0;
0261                         dFlambdadrho = 0;
0262                         dFlambdadtheta = -(1.0 - polmode)*sin(theta)/2.0;
0263                         dFlambdadphi = 0;
0264                         %
0265                         Fdelta = exp(-log(2.0)*(rho - sigmar_rho)^2.0/sigmar_hwhm^2);
0266                         dFdeltadrho = -Fdelta*2.0*log(2.0)*(rho - sigmar_rho)/sigmar_hwhm^2;
0267                         dFdeltadtheta = 0;
0268                         dFdeltadphi = 0;
0269                         %
0270                         ac = 2*sqrt(2*pi)*sigmar_max*sqrt(epsi_rho*epsi_theta/L_rho/xL_perp_fit);
0271                         dacdrho = -0.5*ac*(dL_rhodrho/L_rho + dxL_perpdrho_fit/xL_perp_fit);
0272                         dacdtheta = -0.5*ac*(dL_rhodtheta/L_rho + dxL_perpdtheta_fit/xL_perp_fit);
0273                         dacdphi = -0.5*ac*(dL_rhodphi/L_rho + dxL_perpdphi_fit/xL_perp_fit);                        
0274                         %
0275                         [ht,dht] = halftanh(ac*Flambda*Fdelta*sum(sum(alm.*sin(phase),2),1));
0276                         %
0277                         fluctval.([char(v),'r']) = ht;
0278                         %
0279                         fluctval.(['d',char(v),'rdrho']) = dht*dacdrho*Flambda*Fdelta*sum(sum(alm.*sin(phase),2),1)...
0280                                                          + dht*ac*dFlambdadrho*Fdelta*sum(sum(alm.*sin(phase),2),1)...
0281                                                          + dht*ac*Flambda*dFdeltadrho*sum(sum(alm.*sin(phase),2),1)...
0282                                                          + dht*ac*Flambda*Fdelta*sum(sum(dalmdrho.*sin(phase),2),1)...
0283                                                          + dht*ac*Flambda*Fdelta*sum(sum(alm.*cos(phase).*dphasedrho,2),1);
0284                         %
0285                         fluctval.(['d',char(v),'rdtheta']) = dht*dacdtheta*Flambda*Fdelta*sum(sum(alm.*sin(phase),2),1)...
0286                                                          + dht*ac*dFlambdadtheta*Fdelta*sum(sum(alm.*sin(phase),2),1)...
0287                                                          + dht*ac*Flambda*dFdeltadtheta*sum(sum(alm.*sin(phase),2),1)...
0288                                                          + dht*ac*Flambda*Fdelta*sum(sum(dalmdtheta.*sin(phase),2),1)...
0289                                                          + dht*ac*Flambda*Fdelta*sum(sum(alm.*cos(phase).*dphasedtheta,2),1);
0290                         %
0291                         fluctval.(['d',char(v),'rdphi']) = dht*dacdphi*Flambda*Fdelta*sum(sum(alm.*sin(phase),2),1)...
0292                                                          + dht*ac*dFlambdadphi*Fdelta*sum(sum(alm.*sin(phase),2),1)...
0293                                                          + dht*ac*Flambda*dFdeltadphi*sum(sum(alm.*sin(phase),2),1)...
0294                                                          + dht*ac*Flambda*Fdelta*sum(sum(dalmdphi.*sin(phase),2),1)...
0295                                                          + dht*ac*Flambda*Fdelta*sum(sum(alm.*cos(phase).*dphasedphi,2),1);
0296                         %
0297 
0298                     elseif strcmp(char(field_name(i)),'B_fit'),
0299                         
0300                         
0301                         
0302                         
0303                     end
0304                 elseif fluct_fit.(char(field_name(i)))(ifield).model == 5,%(1) -> 1-D Gaussian drift-wave like model local (cn,cm), absolute epsi values (m)
0305                     if strcmp(char(field_name(i)),'ne_fit'),
0306                         m = [mmin:mmax];
0307                         %
0308                         u_m = pi*epsi_theta*[mmin:mmax]/xL_perp_fit;%here, epsi_perp = epsi_theta
0309                         %
0310                         du_mdrho = -u_m/xL_perp_fit*dxL_perpdrho_fit;
0311                         %
0312                         du_mdtheta = -u_m/xL_perp_fit*dxL_perpdtheta_fit;
0313                         %
0314                         du_mdphi = -u_m/xL_perp_fit*dxL_perpdphi_fit;
0315                         %
0316                         expu_m = exp(-u_m.^2/2);
0317                         am = expu_m;
0318                         %
0319                         damdrho = -expu_m.*u_m.*du_mdrho;
0320                         %
0321                         damdtheta = -expu_m.*u_m.*du_mdtheta;                 
0322                         %
0323                         damdphi = -expu_m.*u_m.*du_mdphi;
0324                         %
0325                         kperp = 2*pi*[mmin:mmax]/xL_perp_fit;
0326                         dkperpdrho = -kperp/xL_perp_fit*dxL_perpdrho_fit;
0327                         dkperpdtheta = 0;
0328                         dkperpdphi = 0;
0329                         %
0330                         phase = (Xicm_fit + Xcn0_fit*(phi-phi0))*kperp + reshape(phaseinit,1,length(mmin:mmax));
0331                         dphasedrho = (dXicmdrho_fit + dXcn0drho_fit*(phi-phi0))*kperp + (Xicm_fit + Xcn0_fit*(phi-phi0))*dkperpdrho;
0332                         dphasedtheta = Xcm_fit*kperp + (Xicm_fit + Xcn0_fit*(phi-phi0))*dkperpdtheta;
0333                         dphasedphi = Xcn0_fit*kperp + (Xicm_fit + Xcn0_fit*(phi-phi0))*dkperpdphi;
0334                         %
0335                         Flambda = (1.0 + polmode + (1.0 - polmode)*cos(theta))/2.0;
0336                         dFlambdadrho = 0;
0337                         dFlambdadtheta = -(1.0 - polmode)*sin(theta)/2.0;
0338                         dFlambdadphi = 0;
0339                         %
0340                         Fdelta = exp(-log(2.0)*(rho - sigmar_rho)^2.0/sigmar_hwhm^2);
0341                         dFdeltadrho = -Fdelta*2.0*log(2.0)*(rho - sigmar_rho)/sigmar_hwhm^2;
0342                         dFdeltadtheta = 0;
0343                         dFdeltadphi = 0;
0344                         %
0345                         ac = 2*sqrt(sqrt(pi))*sigmar_max*sqrt(epsi_theta/xL_perp_fit);
0346                         dacdrho = -0.5*ac*dxL_perpdrho_fit/xL_perp_fit;
0347                         dacdtheta = -0.5*ac*dxL_perpdtheta_fit/xL_perp_fit;
0348                         dacdphi = -0.5*ac*dxL_perpdphi_fit/xL_perp_fit;
0349                         %
0350                         [ht,dht] = halftanh(ac*Flambda*Fdelta*sum(sum(am.*sin(phase),2),1));
0351                         %
0352                         fluctval.([char(v),'r']) = ht;
0353                         %
0354                         fluctval.(['d',char(v),'rdrho']) = dht*dacdrho*Flambda*Fdelta*sum(sum(am.*sin(phase),2),1)...
0355                                                          + dht*ac*dFlambdadrho*Fdelta*sum(sum(am.*sin(phase),2),1)...
0356                                                          + dht*ac*Flambda*dFdeltadrho*sum(sum(am.*sin(phase),2),1)...
0357                                                          + dht*ac*Flambda*Fdelta*sum(sum(damdrho.*sin(phase),2),1)...
0358                                                          + dht*ac*Flambda*Fdelta*sum(sum(am.*cos(phase).*dphasedrho,2),1);
0359                         %
0360                         fluctval.(['d',char(v),'rdtheta']) = dht*dacdtheta*Flambda*Fdelta*sum(sum(am.*sin(phase),2),1)...
0361                                                          + dht*ac*dFlambdadtheta*Fdelta*sum(sum(am.*sin(phase),2),1)...
0362                                                          + dht*ac*Flambda*dFdeltadtheta*sum(sum(am.*sin(phase),2),1)...
0363                                                          + dht*ac*Flambda*Fdelta*sum(sum(damdtheta.*sin(phase),2),1)...
0364                                                          + dht*ac*Flambda*Fdelta*sum(sum(am.*cos(phase).*dphasedtheta,2),1);
0365                         %
0366                         fluctval.(['d',char(v),'rdphi']) = dht*dacdphi*Flambda*Fdelta*sum(sum(am.*sin(phase),2),1)...
0367                                                          + dht*ac*dFlambdadphi*Fdelta*sum(sum(am.*sin(phase),2),1)...
0368                                                          + dht*ac*Flambda*dFdeltadphi*sum(sum(am.*sin(phase),2),1)...
0369                                                          + dht*ac*Flambda*Fdelta*sum(sum(damdphi.*sin(phase),2),1)...
0370                                                          + dht*ac*Flambda*Fdelta*sum(sum(am.*cos(phase).*dphasedphi,2),1);
0371                         %
0372                     elseif strcmp(char(field_name(i)),'B_fit'),
0373                         
0374                         
0375                         
0376                         
0377                     end                    
0378                  elseif fluct_fit.(char(field_name(i)))(ifield).model == 4,%(1) -> 1-D Gaussian drift-wave like model global (rho,curtheta), absolute epsi values (m)
0379                     %
0380                     if strcmp(char(field_name(i)),'ne_fit'),
0381                         
0382                         m = [mmin:mmax];
0383                         %
0384                         u_m = pi*epsi_theta*[mmin:mmax]/xL_perp_fit;%here, epsi_perp = epsi_theta
0385                         %
0386                         du_mdrho = -u_m/xL_perp_fit*dxL_perpdrho_fit;
0387                         %
0388                         du_mdtheta = -u_m/xL_perp_fit*dxL_perpdtheta_fit;
0389                         %
0390                         du_mdphi = -u_m/xL_perp_fit*dxL_perpdphi_fit;
0391                         %
0392                         expu_m = exp(-u_m.^2/2);
0393                         am = expu_m;
0394                         %
0395                         damdrho = -expu_m.*u_m.*du_mdrho;
0396                         %
0397                         damdtheta = -expu_m.*u_m.*du_mdtheta;                 
0398                         %
0399                         damdphi = -expu_m.*u_m.*du_mdphi;
0400                         %
0401                         curtheta = xq_fit*(theta - theta0) - (phi - phi0);
0402                         dcurthetadrho = dxqdrho_fit*(theta - theta0);
0403                         dcurthetadtheta = xq_fit;
0404                         dcurthetadphi = -1;                    
0405                         %
0406                         phase = m*curtheta + reshape(phaseinit,1,length(mmin:mmax));
0407                         dphasedrho = m*dcurthetadrho;
0408                         dphasedtheta = m*dcurthetadtheta;
0409                         dphasedphi = m*dcurthetadphi;
0410                         %
0411                         Flambda = (1.0 + polmode + (1.0 - polmode)*cos(theta))/2.0;
0412                         dFlambdadrho = 0;
0413                         dFlambdadtheta = -(1.0 - polmode)*sin(theta)/2.0;
0414                         dFlambdadphi = 0;
0415                         %
0416                         Fdelta = exp(-log(2.0)*(rho - sigmar_rho)^2.0/sigmar_hwhm^2);
0417                         dFdeltadrho = -Fdelta*2.0*log(2.0)*(rho - sigmar_rho)/sigmar_hwhm^2;
0418                         dFdeltadtheta = 0;
0419                         dFdeltadphi = 0;
0420                         %
0421                         ac = 2*sqrt(sqrt(pi))*sigmar_max*sqrt(epsi_theta/xL_perp_fit);
0422                         dacdrho = -0.5*ac*dxL_perpdrho_fit/xL_perp_fit;
0423                         dacdtheta = -0.5*ac*dxL_perpdtheta_fit/xL_perp_fit;
0424                         dacdphi = -0.5*ac*dxL_perpdphi_fit/xL_perp_fit;
0425                         %
0426                         [ht,dht] = halftanh(ac*Flambda*Fdelta*sum(sum(am.*sin(phase),2),1));
0427                         %
0428                         fluctval.([char(v),'r']) = ht;
0429                         %
0430                         fluctval.(['d',char(v),'rdrho']) = dht*dacdrho*Flambda*Fdelta*sum(sum(am.*sin(phase),2),1)...
0431                                                          + dht*ac*dFlambdadrho*Fdelta*sum(sum(am.*sin(phase),2),1)...
0432                                                          + dht*ac*Flambda*dFdeltadrho*sum(sum(am.*sin(phase),2),1)...
0433                                                          + dht*ac*Flambda*Fdelta*sum(sum(damdrho.*sin(phase),2),1)...
0434                                                          + dht*ac*Flambda*Fdelta*sum(sum(am.*cos(phase).*dphasedrho,2),1);
0435                         %
0436                         fluctval.(['d',char(v),'rdtheta']) = dht*dacdtheta*Flambda*Fdelta*sum(sum(am.*sin(phase),2),1)...
0437                                                          + dht*ac*dFlambdadtheta*Fdelta*sum(sum(am.*sin(phase),2),1)...
0438                                                          + dht*ac*Flambda*dFdeltadtheta*sum(sum(am.*sin(phase),2),1)...
0439                                                          + dht*ac*Flambda*Fdelta*sum(sum(damdtheta.*sin(phase),2),1)...
0440                                                          + dht*ac*Flambda*Fdelta*sum(sum(am.*cos(phase).*dphasedtheta,2),1);
0441                         %
0442                         fluctval.(['d',char(v),'rdphi']) = dht*dacdphi*Flambda*Fdelta*sum(sum(am.*sin(phase),2),1)...
0443                                                          + dht*ac*dFlambdadphi*Fdelta*sum(sum(am.*sin(phase),2),1)...
0444                                                          + dht*ac*Flambda*dFdeltadphi*sum(sum(am.*sin(phase),2),1)...
0445                                                          + dht*ac*Flambda*Fdelta*sum(sum(damdphi.*sin(phase),2),1)...
0446                                                          + dht*ac*Flambda*Fdelta*sum(sum(am.*cos(phase).*dphasedphi,2),1);
0447                         %
0448                     elseif strcmp(char(field_name(i)),'B_fit'),
0449                         
0450                         
0451                         
0452                         
0453                     end                  
0454 
0455 
0456                 elseif fluct_fit.(char(field_name(i)))(ifield).model == 3,%(1) -> 2-D Gaussian drift-wave like model global (rho,curtheta), absolute epsi values (m)
0457                     %
0458                     if strcmp(char(field_name(i)),'ne_fit'),
0459                         %
0460                         [l,m] = build_tensor_yp(lmin:lmax,mmin:mmax);
0461                         %
0462                         u_l = pi*epsi_rho*[lmin:lmax]/L_rho;
0463                         u_m = pi*epsi_theta*[mmin:mmax]/xL_perp_fit;%here, epsi_perp = epsi_theta
0464                         %
0465                         du_ldrho = -u_l/L_rho*dL_rhodrho;
0466                         du_mdrho = -u_m/xL_perp_fit*dxL_perpdrho_fit;
0467                         %
0468                         du_ldtheta = -u_l/L_rho*dL_rhodtheta;
0469                         du_mdtheta = -u_m/xL_perp_fit*dxL_perpdtheta_fit;
0470                         %
0471                         du_ldphi = -u_l/L_rho*dL_rhodphi;
0472                         du_mdphi = -u_m/xL_perp_fit*dxL_perpdphi_fit;
0473                         %
0474                         expu_l = exp(-u_l.^2/2);
0475                         expu_m = exp(-u_m.^2/2);
0476                         %
0477                         [al,am] = build_tensor_yp(expu_l,expu_m);
0478                         alm = al.*am;
0479                         %
0480                         dexpu_ldrho = -expu_l.*u_l.*du_ldrho;
0481                         dexpu_mdrho = -expu_m.*u_m.*du_mdrho;
0482                         %
0483                         [daldrho,damdrho] = build_tensor_yp(dexpu_ldrho,dexpu_mdrho);
0484                         dalmdrho = daldrho.*am + al.*damdrho;
0485                         %
0486                         dexpu_ldtheta = -expu_l.*u_l.*du_ldtheta;
0487                         dexpu_mdtheta = -expu_m.*u_m.*du_mdtheta;
0488                         %
0489                         [daldtheta,damdtheta] = build_tensor_yp(dexpu_ldtheta,dexpu_mdtheta);
0490                         dalmdtheta = daldtheta.*am + al.*damdtheta;                        
0491                         %
0492                         dexpu_ldphi = -expu_l.*u_l.*du_ldphi;
0493                         dexpu_mdphi = -expu_m.*u_m.*du_mdphi;
0494                         %
0495                         [daldphi,damdphi] = build_tensor_yp(dexpu_ldphi,dexpu_mdphi);
0496                         dalmdphi = daldphi.*am + al.*damdphi;
0497                         %
0498                         curtheta = xq_fit*theta - phi;
0499                         dcurthetadrho = dxqdrho_fit*theta;
0500                         dcurthetadtheta = xq_fit;
0501                         dcurthetadphi = -1;                    
0502                         %
0503                         phase = 2*pi*l*rho + m*curtheta + reshape(phaseinit,length(lmin:lmax),length(mmin:mmax));
0504                         dphasedrho = 2*pi*l + m*dcurthetadrho;
0505                         dphasedtheta = m*dcurthetadtheta;
0506                         dphasedphi = m*dcurthetadphi;
0507                         %
0508                         Flambda = (1.0 + polmode + (1.0 - polmode)*cos(theta))/2.0;
0509                         dFlambdadrho = 0;
0510                         dFlambdadtheta = -(1.0 - polmode)*sin(theta)/2.0;
0511                         dFlambdadphi = 0;
0512                         %
0513                         Fdelta = exp(-log(2.0)*(rho - sigmar_rho)^2.0/sigmar_hwhm^2);
0514                         dFdeltadrho = -Fdelta*2.0*log(2.0)*(rho - sigmar_rho)/sigmar_hwhm^2;
0515                         dFdeltadtheta = 0;
0516                         dFdeltadphi = 0;
0517                         %
0518                         ac = 2*sqrt(2*pi)*sigmar_max*sqrt(epsi_rho*epsi_theta/L_rho/xL_perp_fit);
0519                         dacdrho = -0.5*ac*(dL_rhodrho/L_rho + dxL_perpdrho_fit/xL_perp_fit);
0520                         dacdtheta = -0.5*ac*(dL_rhodtheta/L_rho + dxL_perpdtheta_fit/xL_perp_fit);
0521                         dacdphi = -0.5*ac*(dL_rhodphi/L_rho + dxL_perpdphi_fit/xL_perp_fit);
0522                         %
0523                         [ht,dht] = halftanh(ac*Flambda*Fdelta*sum(sum(alm.*sin(phase),2),1));
0524                         %
0525                         fluctval.([char(v),'r']) = ht;
0526                         %
0527                         fluctval.(['d',char(v),'rdrho']) = dht*dacdrho*Flambda*Fdelta*sum(sum(alm.*sin(phase),2),1)...
0528                                                          + dht*ac*dFlambdadrho*Fdelta*sum(sum(alm.*sin(phase),2),1)...
0529                                                          + dht*ac*Flambda*dFdeltadrho*sum(sum(alm.*sin(phase),2),1)...
0530                                                          + dht*ac*Flambda*Fdelta*sum(sum(dalmdrho.*sin(phase),2),1)...
0531                                                          + dht*ac*Flambda*Fdelta*sum(sum(alm.*cos(phase).*dphasedrho,2),1);
0532                         %
0533                         fluctval.(['d',char(v),'rdtheta']) = dht*dacdtheta*Flambda*Fdelta*sum(sum(alm.*sin(phase),2),1)...
0534                                                          + dht*ac*dFlambdadtheta*Fdelta*sum(sum(alm.*sin(phase),2),1)...
0535                                                          + dht*ac*Flambda*dFdeltadtheta*sum(sum(alm.*sin(phase),2),1)...
0536                                                          + dht*ac*Flambda*Fdelta*sum(sum(dalmdtheta.*sin(phase),2),1)...
0537                                                          + dht*ac*Flambda*Fdelta*sum(sum(alm.*cos(phase).*dphasedtheta,2),1);
0538                         %
0539                         fluctval.(['d',char(v),'rdphi']) = dht*dacdphi*Flambda*Fdelta*sum(sum(alm.*sin(phase),2),1)...
0540                                                          + dht*ac*dFlambdadphi*Fdelta*sum(sum(alm.*sin(phase),2),1)...
0541                                                          + dht*ac*Flambda*dFdeltadphi*sum(sum(alm.*sin(phase),2),1)...
0542                                                          + dht*ac*Flambda*Fdelta*sum(sum(dalmdphi.*sin(phase),2),1)...
0543                                                          + dht*ac*Flambda*Fdelta*sum(sum(alm.*cos(phase).*dphasedphi,2),1);
0544                         %
0545                     elseif strcmp(char(field_name(i)),'B_fit'),
0546                         
0547                         
0548                         
0549                         
0550                     end                  
0551                 elseif fluct_fit.(char(field_name(i)))(ifield).model == 2,%(1) -> 3-D Gaussian model (rho,theta,phi), absolute epsi values (benchmark of C3PO)
0552                     %
0553                     if strcmp(char(field_name(i)),'ne_fit'),
0554                         %
0555                         [l,m,n] = build_tensor_yp(lmin:lmax,mmin:mmax,nmin:nmax);
0556                         %
0557                         u_l = pi*epsi_rho*[lmin:lmax]/L_rho;
0558                         u_m = pi*epsi_theta*[mmin:mmax]/xL_theta_fit;
0559                         u_n = pi*epsi_phi*[nmin:nmax]/XL_phi_fit;
0560                         %
0561                         du_ldrho = -u_l/L_rho*dL_rhodrho;
0562                         du_mdrho = -u_m/xL_theta_fit*dxL_thetadrho_fit;
0563                         du_ndrho = -u_n/XL_phi_fit*dXL_phidrho_fit;
0564                         %
0565                         du_ldtheta = -u_l/L_rho*dL_rhodtheta;
0566                         du_mdtheta = -u_m/xL_theta_fit*dxL_thetadtheta_fit;
0567                         du_ndtheta = -u_n/XL_phi_fit*dXL_phidtheta_fit;
0568                         %
0569                         du_ldphi = -u_l/L_rho*dL_rhodphi;
0570                         du_mdphi = -u_m/xL_theta_fit*dxL_thetadphi_fit;
0571                         du_ndphi = -u_n/XL_phi_fit*dXL_phidphi_fit;
0572                         %
0573                         expu_l = exp(-u_l.^2/2);
0574                         expu_m = exp(-u_m.^2/2);
0575                         expu_n = exp(-u_n.^2/2);
0576                         %
0577                         [al,am,an] = build_tensor_yp(expu_l,expu_m,expu_n);
0578                         almn = al.*am.*an;
0579                         %
0580                         dexpu_ldrho = -expu_l.*u_l.*du_ldrho;
0581                         dexpu_mdrho = -expu_m.*u_m.*du_mdrho;
0582                         dexpu_ndrho = -expu_n.*u_n.*du_ndrho;
0583                         %
0584                         [daldrho,damdrho,dandrho] = build_tensor_yp(dexpu_ldrho,dexpu_mdrho,dexpu_ndrho);
0585                         dalmndrho = daldrho.*am.*an + al.*damdrho.*an + al.*am.*dandrho;
0586                         %
0587                         dexpu_ldtheta = -expu_l.*u_l.*du_ldtheta;
0588                         dexpu_mdtheta = -expu_m.*u_m.*du_mdtheta;
0589                         dexpu_ndtheta = -expu_n.*u_n.*du_ndtheta; 
0590                         %
0591                         [daldtheta,damdtheta,dandtheta] = build_tensor_yp(dexpu_ldtheta,dexpu_mdtheta,dexpu_ndtheta);
0592                         dalmndtheta = daldtheta.*am.*an + al.*damdtheta.*an + al.*am.*dandtheta;                        
0593                         %
0594                         dexpu_ldphi = -expu_l.*u_l.*du_ldphi;
0595                         dexpu_mdphi = -expu_m.*u_m.*du_mdphi;
0596                         dexpu_ndphi = -expu_n.*u_n.*du_ndphi;
0597                         %
0598                         [daldphi,damdphi,dandphi] = build_tensor_yp(dexpu_ldphi,dexpu_mdphi,dexpu_ndphi);
0599                         dalmndphi = daldphi.*am.*an + al.*damdphi.*an + al.*am.*dandphi;
0600                         %
0601                         phase = 2*pi*l*rho + m*theta + n*phi + reshape(phaseinit,length(lmin:lmax),length(mmin:mmax),length(nmin:nmax));
0602                         dphasedrho = 2*pi*l;
0603                         dphasedtheta = m;
0604                         dphasedphi = n;
0605                         %
0606                         Flambda = (1.0 + polmode + (1.0 - polmode)*cos(theta))/2.0;
0607                         dFlambdadrho = 0;
0608                         dFlambdadtheta = -(1.0 - polmode)*sin(theta)/2.0;
0609                         dFlambdadphi = 0;
0610                         %
0611                         Fdelta = exp(-log(2.0)*(rho - sigmar_rho)^2.0/sigmar_hwhm^2);
0612                         dFdeltadrho = -Fdelta*2.0*log(2.0)*(rho - sigmar_rho)/sigmar_hwhm^2;
0613                         dFdeltadtheta = 0;
0614                         dFdeltadphi = 0;
0615                         %
0616                         ac = 4*pi^0.75*sigmar_max*sqrt(epsi_rho*epsi_theta*epsi_phi/L_rho/xL_theta_fit/XL_phi_fit);
0617                         dacdrho = -0.5*ac*(dL_rhodrho/L_rho + dxL_thetadrho_fit/xL_theta_fit + dXL_phidrho_fit/XL_phi_fit);
0618                         dacdtheta = -0.5*ac*(dL_rhodtheta/L_rho + dxL_thetadtheta_fit/xL_theta_fit + dXL_phidtheta_fit/XL_phi_fit);
0619                         dacdphi = -0.5*ac*(dL_rhodphi/L_rho + dxL_thetadphi_fit/xL_theta_fit +  dXL_phidphi_fit/XL_phi_fit);
0620                         %
0621                         [ht,dht] = halftanh(ac*Flambda*Fdelta*sum(sum(sum(almn.*sin(phase),3),2),1));
0622                         %
0623                         fluctval.([char(v),'r']) = ht;
0624                         %
0625                         fluctval.(['d',char(v),'rdrho']) = dht*dacdrho*Flambda*Fdelta*sum(sum(sum(almn.*sin(phase),3),2),1)...
0626                                                          + dht*ac*dFlambdadrho*Fdelta*sum(sum(sum(almn.*sin(phase),3),2),1)...
0627                                                          + dht*ac*Flambda*dFdeltadrho*sum(sum(sum(almn.*sin(phase),3),2),1)...
0628                                                          + dht*ac*Flambda*Fdelta*sum(sum(sum(dalmndrho.*sin(phase),3),2),1)...
0629                                                          + dht*ac*Flambda*Fdelta*sum(sum(sum(almn.*cos(phase).*dphasedrho,3),2),1);
0630                         %
0631                         fluctval.(['d',char(v),'rdtheta']) = dht*dacdtheta*Flambda*Fdelta*sum(sum(sum(almn.*sin(phase),3),2),1)...
0632                                                          + dht*ac*dFlambdadtheta*Fdelta*sum(sum(sum(almn.*sin(phase),3),2),1)...
0633                                                          + dht*ac*Flambda*dFdeltadtheta*sum(sum(sum(almn.*sin(phase),3),2),1)...
0634                                                          + dht*ac*Flambda*Fdelta*sum(sum(sum(dalmndtheta.*sin(phase),3),2),1)...
0635                                                          + dht*ac*Flambda*Fdelta*sum(sum(sum(almn.*cos(phase).*dphasedtheta,3),2),1);
0636                         %
0637                         fluctval.(['d',char(v),'rdphi']) = dht*dacdphi*Flambda*Fdelta*sum(sum(sum(almn.*sin(phase),3),2),1)...
0638                                                          + dht*ac*dFlambdadphi*Fdelta*sum(sum(sum(almn.*sin(phase),3),2),1)...
0639                                                          + dht*ac*Flambda*dFdeltadphi*sum(sum(sum(almn.*sin(phase),3),2),1)...
0640                                                          + dht*ac*Flambda*Fdelta*sum(sum(sum(dalmndphi.*sin(phase),3),2),1)...
0641                                                          + dht*ac*Flambda*Fdelta*sum(sum(sum(almn.*cos(phase).*dphasedphi,3),2),1);
0642                         %
0643                     elseif strcmp(char(field_name(i)),'B_fit'),
0644                         
0645                         
0646                         
0647                         
0648                     end
0649                 elseif fluct_fit.(char(field_name(i)))(ifield).model == 1,%(1) -> 3-D Gaussian model (rho,theta,phi), relative epsi values (benchmark of C3PO)
0650                     %
0651                     %dkperp = kperp(2) - kperp(1);%Uniform kperp grid assumed
0652                     %
0653                     %keyboard
0654                     if strcmp(char(field_name(i)),'ne_fit'),
0655 
0656                         [l,m,n] = build_tensor_yp(lmin:lmax,mmin:mmax,nmin:nmax);
0657                         %
0658                         u_l = pi*epsi_rho*[lmin:lmax];
0659                         u_m = pi*epsi_theta*[mmin:mmax];
0660                         u_n = pi*epsi_phi*[nmin:nmax];
0661                         %
0662                         expu_l = exp(-u_l.^2/2);
0663                         expu_m = exp(-u_m.^2/2);
0664                         expu_n = exp(-u_n.^2/2);
0665                         %
0666                         [al,am,an] = build_tensor_yp(expu_l,expu_m,expu_n);
0667                         almn = al.*am.*an;                        
0668                         %
0669                         phase = 2*pi*l*rho + m*theta + n*phi + reshape(phaseinit,length(lmin:lmax),length(mmin:mmax),length(nmin:nmax));
0670                         dphasedrho = 2*pi*l;
0671                         dphasedtheta = m;
0672                         dphasedphi = n;
0673                         %
0674                         Flambda = (1.0 + polmode + (1.0 - polmode)*cos(theta))/2.0;
0675                         dFlambdadtheta = -(1.0 - polmode)*sin(theta)/2.0;
0676                         %
0677                         Fdelta = exp(-log(2.0)*(rho - sigmar_rho)^2.0/sigmar_hwhm^2);
0678                         dFdeltadrho = -Fdelta*2.0*log(2.0)*(rho - sigmar_rho)/sigmar_hwhm^2;
0679                         %
0680                         ac = 4*pi^0.75*sigmar_max*sqrt(epsi_rho*epsi_theta*epsi_phi);
0681                         %
0682                         [ht,dht] = halftanh(ac*Flambda*Fdelta*sum(sum(sum(almn.*sin(phase),3),2),1));
0683                         %
0684                         fluctval.([char(v),'r']) = ht;
0685                         %
0686                         fluctval.(['d',char(v),'rdrho']) = dht*ac*Flambda*Fdelta*sum(sum(sum(almn.*cos(phase).*dphasedrho,3),2),1)...
0687                                                          + dht*ac*Flambda*dFdeltadrho*sum(sum(sum(almn.*sin(phase),3),2),1);
0688                         %
0689                         fluctval.(['d',char(v),'rdtheta']) = dht*ac*Flambda*Fdelta*sum(sum(sum(almn.*cos(phase).*dphasedtheta,3),2),1)...
0690                                                            + dht*ac*dFlambdadtheta*Fdelta*sum(sum(sum(almn.*sin(phase),3),2),1);
0691                         %
0692                         fluctval.(['d',char(v),'rdphi']) = dht*ac*Flambda*Fdelta*sum(sum(sum(almn.*cos(phase).*dphasedphi,3),2),1);
0693                                              
0694                         %
0695                         %for ikperp = 1:nkperp,
0696                             %
0697                        %     phase = gphase*kperp(ikperp) + phaseinit(ikperp);
0698                         %    dphasedrho = dgphasedrho.*kperp(ikperp);
0699                        %     dphasedtheta = dgphasedtheta.*kperp(ikperp);
0700                        %     dphasedphi = dgphasedphi.*kperp(ikperp);
0701                        %     %
0702                        %     df(ikperp,:) = w_fit{1}.*cos(phase).*sqrt(2.*dkperp./w_fit{2}).*exp(-(kperp(ikperp) - w_fit{3}).^2./w_fit{2}./w_fit{2}./4)/(2*pi)^(0.25);
0703                             %
0704                        %     ddfdrho(ikperp,:) = dwdrho_fit{1}.*cos(phase).*sqrt(2.*dkperp./w_fit{2}).*exp(-(kperp(ikperp) - w_fit{3}).^2./w_fit{2}./w_fit{2}./4)/(2*pi)^(0.25)...
0705                        %                         - w_fit{1}.*dphasedrho.*sin(phase).*sqrt(2.*dkperp./w_fit{2}).*exp(-(kperp(ikperp) - w_fit{3}).^2./w_fit{2}./w_fit{2}./4)/(2*pi)^(0.25)...
0706                        %                         - 0.5*w_fit{1}.*cos(phase).*dwdrho_fit{2}.*sqrt(2.*dkperp./w_fit{2})./w_fit{2}.*exp(-(kperp(ikperp) - w_fit{3}).^2./w_fit{2}./w_fit{2}./4)/(2*pi)^(0.25)...
0707                        %                         + w_fit{1}.*cos(phase).*sqrt(2.*dkperp./w_fit{2}).*exp(-(kperp(ikperp) - w_fit{3}).^2./w_fit{2}./w_fit{2}./4)/(2*pi)^(0.25).*(dwdrho_fit{3}.*w_fit{2} + (kperp(ikperp) - w_fit{3}).*dwdrho_fit{2}).*(kperp(ikperp) - w_fit{3})./w_fit{2}.^3./2;
0708                             %
0709                        %     ddfdtheta(ikperp,:) = dwdtheta_fit{1}.*cos(phase).*sqrt(2.*dkperp./w_fit{2}).*exp(-(kperp(ikperp) - w_fit{3}).^2./w_fit{2}./w_fit{2}./4)/(2*pi)^(0.25)...
0710                        %                         - w_fit{1}.*dphasedtheta.*sin(phase).*sqrt(2.*dkperp./w_fit{2}).*exp(-(kperp(ikperp) - w_fit{3}).^2./w_fit{2}./w_fit{2}./4)/(2*pi)^(0.25)...
0711                        %                         - 0.5*w_fit{1}.*cos(phase).*dwdtheta_fit{2}.*sqrt(2.*dkperp./w_fit{2})./w_fit{2}.*exp(-(kperp(ikperp) - w_fit{3}).^2./w_fit{2}./w_fit{2}./4)/(2*pi)^(0.25)...
0712                         %                        + w_fit{1}.*cos(phase).*sqrt(2.*dkperp./w_fit{2}).*exp(-(kperp(ikperp) - w_fit{3}).^2./w_fit{2}./w_fit{2}./4)/(2*pi)^(0.25).*(dwdtheta_fit{3}.*w_fit{2} + (kperp(ikperp) - w_fit{3}).*dwdtheta_fit{2}).*(kperp(ikperp) - w_fit{3})./w_fit{2}.^3./2;
0713                             %
0714                         %    ddfdphi(ikperp,:) = -w_fit{1}.*dphasedphi.*sin(phase).*sqrt(2.*dkperp./w_fit{2}).*exp(-(kperp(ikperp) - w_fit{3}).^2./w_fit{2}./w_fit{2}./4)/(2*pi)^(0.25);
0715                         %end
0716                         %
0717                         %fluctval.([char(v),'r']) = fluctval.([char(v),'r']) + sum(df,1);
0718                         %fluctval.(['d',char(v),'rdrho']) = fluctval.(['d',char(v),'rdrho']) + sum(ddfdrho,1);
0719                         %fluctval.(['d',char(v),'rdtheta']) = fluctval.(['d',char(v),'rdtheta']) + sum(ddfdtheta,1);
0720                         %fluctval.(['d',char(v),'rdphi']) = fluctval.(['d',char(v),'rdphi']) + sum(ddfdphi,1);
0721                         %
0722 
0723                         %
0724                     elseif strcmp(char(field_name(i)),'B_fit')
0725                         dfx = zeros(nkperp,length(theta));
0726                         ddfxdrho = zeros(nkperp,length(theta));
0727                         ddfxdtheta = zeros(nkperp,length(theta));
0728                         dfxdphi = zeros(nkperp,length(theta));
0729                         %
0730                         dfy = zeros(nkperp,length(theta));
0731                         ddfydrho = zeros(nkperp,length(theta));
0732                         ddfydtheta = zeros(nkperp,length(theta));
0733                         dfydphi = zeros(nkperp,length(theta));
0734                         %
0735                         dfz = zeros(nkperp,length(theta));
0736                         ddfzdrho = zeros(nkperp,length(theta));
0737                         ddfzdtheta = zeros(nkperp,length(theta));
0738                         dfzdphi = zeros(nkperp,length(theta));
0739                         %
0740                         %for ikperp = 1:nkperp,
0741                             %
0742                         %    phase = gphase*kperp(ikperp) + phaseinit(ikperp);
0743                          %   dphasedrho = dgphasedrho.*kperp(ikperp);
0744                         %    dphasedtheta = dgphasedtheta.*kperp(ikperp);
0745                         %    dphasedphi = dgphasedphi.*kperp(ikperp);
0746                             %
0747                         %    dfx(ikperp,:) = w_fit{1}.*cos(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25);
0748                             %
0749                          %   ddfxdrho(ikperp,:) = dwdrho_fit{1}.*cos(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0750                           %                      - w_fit{1}.*dphasedrho.*sin(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0751                          %                       - 0.5*w_fit{1}.*cos(phase).*dwdrho_fit{2}.*sqrt(2.*dkperp./w_fit{4})./w_fit{4}.*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0752                          %                       + w_fit{1}.*cos(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25).*(dwdrho_fit{3}.*w_fit{4} + (kperp(ikperp) - w_fit{5}).*dwdrho_fit{2}).*(kperp(ikperp) - w_fit{5})./w_fit{4}.^3./2;
0753                             %
0754                          %   ddfxdtheta(ikperp,:) = dwdtheta_fit{1}.*cos(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0755                           %                      - w_fit{1}.*dphasedtheta.*sin(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0756                           %                      - 0.5*w_fit{1}.*cos(phase).*dwdtheta_fit{2}.*sqrt(2.*dkperp./w_fit{4})./w_fit{4}.*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0757                            %                     + w_fit{1}.*cos(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25).*(dwdtheta_fit{3}.*w_fit{4} + (kperp(ikperp) - w_fit{5}).*dwdtheta_fit{2}).*(kperp(ikperp) - w_fit{5})./w_fit{4}.^3./2;
0758                             %
0759                           %  ddfxdphi(ikperp,:) = -w_fit{1}.*dphasedphi.*sin(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25);
0760                             %
0761                           %  dfy(ikperp,:) = w_fit{2}.*cos(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25);
0762                             %
0763                           %  ddfydrho(ikperp,:) = dwdrho_fit{1}.*cos(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0764                            %                     - w_fit{2}.*dphasedrho.*sin(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0765                           %                      - 0.5*w_fit{2}.*cos(phase).*dwdrho_fit{2}.*sqrt(2.*dkperp./w_fit{4})./w_fit{4}.*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0766                            %                     + w_fit{2}.*cos(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25).*(dwdrho_fit{3}.*w_fit{4} + (kperp(ikperp) - w_fit{5}).*dwdrho_fit{2}).*(kperp(ikperp) - w_fit{5})./w_fit{4}.^3./2;
0767                             %
0768                            % ddfydtheta(ikperp,:) = dwdtheta_fit{1}.*cos(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0769                            %                     - w_fit{2}.*dphasedtheta.*sin(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0770                            %                     - 0.5*w_fit{2}.*cos(phase).*dwdtheta_fit{2}.*sqrt(2.*dkperp./w_fit{4})./w_fit{4}.*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0771                            %                     + w_fit{2}.*cos(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25).*(dwdtheta_fit{3}.*w_fit{4} + (kperp(ikperp) - w_fit{5}).*dwdtheta_fit{2}).*(kperp(ikperp) - w_fit{5})./w_fit{4}.^3./2;
0772                             %
0773                            % ddfydphi(ikperp,:) = -w_fit{2}.*dphasedphi.*sin(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25);
0774                             %
0775                             %
0776                            % dfz(ikperp,:) = w_fit{3}.*cos(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25);
0777                             %
0778                            % ddfzdrho(ikperp,:) = dwdrho_fit{1}.*cos(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0779                            %                     - w_fit{3}.*dphasedrho.*sin(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0780                            %                     - 0.5*w_fit{3}.*cos(phase).*dwdrho_fit{2}.*sqrt(2.*dkperp./w_fit{4})./w_fit{4}.*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0781                            %                     + w_fit{3}.*cos(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25).*(dwdrho_fit{3}.*w_fit{4} + (kperp(ikperp) - w_fit{5}).*dwdrho_fit{2}).*(kperp(ikperp) - w_fit{5})./w_fit{4}.^3./2;
0782                             %
0783                            % ddfzdtheta(ikperp,:) = dwdtheta_fit{1}.*cos(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0784                            %                     - w_fit{3}.*dphasedtheta.*sin(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0785                            %                     - 0.5*w_fit{3}.*cos(phase).*dwdtheta_fit{2}.*sqrt(2.*dkperp./w_fit{4})./w_fit{4}.*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25)...
0786                            %                     + w_fit{3}.*cos(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25).*(dwdtheta_fit{3}.*w_fit{4} + (kperp(ikperp) - w_fit{5}).*dwdtheta_fit{2}).*(kperp(ikperp) - w_fit{5})./w_fit{4}.^3./2;
0787                             %
0788                            % ddfzdphi(ikperp,:) = -w_fit{3}.*dphasedphi.*sin(phase).*sqrt(2.*dkperp./w_fit{4}).*exp(-(kperp(ikperp) - w_fit{5}).^2./w_fit{4}./w_fit{4}./4)/(2*pi)^(0.25);
0789                             %
0790                         %end
0791                         %
0792                         %fluctval.([char(v),'xr']) = fluctval.([char(v),'xr']) + sum(dfx,1);
0793                         %fluctval.(['d',char(v),'xrdrho']) = fluctval.(['d',char(v),'xrdrho']) + sum(ddfxdrho,1);
0794                         %fluctval.(['d',char(v),'xrdtheta']) = fluctval.(['d',char(v),'xrdtheta']) + sum(ddfxdtheta,1);
0795                         %fluctval.(['d',char(v),'xrdphi']) = fluctval.(['d',char(v),'xrdphi']) + sum(ddfxdphi,1);
0796                         %
0797                         %fluctval.([char(v),'yr']) = fluctval.([char(v),'yr']) + sum(dfy,1);
0798                         %fluctval.(['d',char(v),'yrdrho']) = fluctval.(['d',char(v),'yrdrho']) + sum(ddfydrho,1);
0799                         %fluctval.(['d',char(v),'yrdtheta']) = fluctval.(['d',char(v),'yrdtheta']) + sum(ddfydtheta,1);
0800                         %fluctval.(['d',char(v),'yrdphi']) = fluctval.(['d',char(v),'yrdphi']) + sum(ddfydphi,1);
0801                         %
0802                         %fluctval.([char(v),'zr']) = fluctval.([char(v),'zr']) + sum(dfz,1);
0803                         %fluctval.(['d',char(v),'zrdrho']) = fluctval.(['d',char(v),'zrdrho']) + sum(ddfzdrho,1);
0804                         %fluctval.(['d',char(v),'zrdtheta']) = fluctval.(['d',char(v),'zrdtheta']) + sum(ddfzdtheta,1);
0805                         %fluctval.(['d',char(v),'zrdphi']) = fluctval.(['d',char(v),'zrdphi']) + sum(ddfzdphi,1);
0806                         %
0807                         fluctval.([char(v),'xr']) = fluctval.([char(v),'xr']) + sum(dfx,1);
0808                         fluctval.(['d',char(v),'xrdrho']) = fluctval.(['d',char(v),'xrdrho']) + sum(ddfxdrho,1);
0809                         fluctval.(['d',char(v),'xrdtheta']) = fluctval.(['d',char(v),'xrdtheta']) + sum(ddfxdtheta,1);
0810                         fluctval.(['d',char(v),'xrdphi']) = fluctval.(['d',char(v),'xrdphi']) + sum(ddfxdphi,1);
0811                         %
0812                         fluctval.([char(v),'yr']) = fluctval.([char(v),'yr']) + sum(dfy,1);
0813                         fluctval.(['d',char(v),'yrdrho']) = fluctval.(['d',char(v),'yrdrho']) + sum(ddfydrho,1);
0814                         fluctval.(['d',char(v),'yrdtheta']) = fluctval.(['d',char(v),'yrdtheta']) + sum(ddfydtheta,1);
0815                         fluctval.(['d',char(v),'yrdphi']) = fluctval.(['d',char(v),'yrdphi']) + sum(ddfydphi,1);
0816                         %
0817                         fluctval.([char(v),'zr']) = fluctval.([char(v),'zr']) + sum(dfz,1);
0818                         fluctval.(['d',char(v),'zrdrho']) = fluctval.(['d',char(v),'zrdrho']) + sum(ddfzdrho,1);
0819                         fluctval.(['d',char(v),'zrdtheta']) = fluctval.(['d',char(v),'zrdtheta']) + sum(ddfzdtheta,1);
0820                         fluctval.(['d',char(v),'zrdphi']) = fluctval.(['d',char(v),'zrdphi']) + sum(ddfzdphi,1);
0821                         %
0822                     end
0823                 elseif fluct_fit.(char(field_name(i)))(ifield).model == 0,%Magnetic ripple model
0824                     
0825                     % Magnetic field
0826                                     
0827                     dfx = zeros(1,length(theta));
0828                     ddfxdrho = zeros(1,length(theta));
0829                     ddfxdtheta = zeros(1,length(theta));
0830                     dfxdphi = zeros(1,length(theta));
0831                     %
0832                     dfy = zeros(1,length(theta));
0833                     ddfydrho = zeros(1,length(theta));
0834                     ddfydtheta = zeros(1,length(theta));
0835                     dfydphi = zeros(1,length(theta));
0836                     %
0837                     dfz = zeros(1,length(theta));
0838                     ddfzdrho = zeros(1,length(theta));
0839                     ddfzdtheta = zeros(1,length(theta));
0840                     dfzdphi = zeros(1,length(theta));
0841                     %
0842                     dfx(1) = w_fit{1}.*sin(nmin*phi);
0843                     ddfxdrho(1) = dwdrho_fit{1}.*sin(nmin*phi);                  
0844                     ddfxdtheta(1) = dwdtheta_fit{1}.*sin(nmin*phi);                  
0845                     ddfxdphi(1) = w_fit{1}.*cos(nmin*phi).*nmin; 
0846                     %
0847                     dfy(1) = w_fit{2}.*sin(nmin*phi);
0848                     ddfydrho(1) = dwdrho_fit{2}.*sin(nmin*phi);                  
0849                     ddfydtheta(1) = dwdtheta_fit{2}.*sin(nmin*phi);                  
0850                     ddfydphi(1) = w_fit{2}.*cos(nmin*phi).*nmin; 
0851                     %
0852                     dfz(1) = w_fit{3}.*cos(nmin*phi);
0853                     ddfzdrho(1) = dwdrho_fit{3}.*cos(nmin*phi);                  
0854                     ddfzdtheta(1) = dwdtheta_fit{3}.*cos(nmin*phi);                  
0855                     ddfzdphi(1) = -w_fit{3}.*sin(nmin*phi).*nmin; 
0856                     %
0857                     fluctval.([char(v),'xr']) = fluctval.([char(v),'xr']) + dfx;
0858                     fluctval.(['d',char(v),'xrdrho']) = fluctval.(['d',char(v),'xrdrho']) + ddfxdrho;
0859                     fluctval.(['d',char(v),'xrdtheta']) = fluctval.(['d',char(v),'xrdtheta']) + ddfxdtheta;
0860                     fluctval.(['d',char(v),'xrdphi']) = fluctval.(['d',char(v),'xrdphi']) + ddfxdphi;
0861                     %
0862                     fluctval.([char(v),'yr']) = fluctval.([char(v),'yr']) + dfy;
0863                     fluctval.(['d',char(v),'yrdrho']) = fluctval.(['d',char(v),'yrdrho']) + ddfydrho;
0864                     fluctval.(['d',char(v),'yrdtheta']) = fluctval.(['d',char(v),'yrdtheta']) + ddfydtheta;
0865                     fluctval.(['d',char(v),'yrdphi']) = fluctval.(['d',char(v),'yrdphi']) + ddfydphi;
0866                     %
0867                     fluctval.([char(v),'zr']) = fluctval.([char(v),'zr']) + dfz;
0868                     fluctval.(['d',char(v),'zrdrho']) = fluctval.(['d',char(v),'zrdrho']) + ddfzdrho;
0869                     fluctval.(['d',char(v),'zrdtheta']) = fluctval.(['d',char(v),'zrdtheta']) + ddfzdtheta;
0870                     fluctval.(['d',char(v),'zrdphi']) = fluctval.(['d',char(v),'zrdphi']) + ddfzdphi;
0871                     %
0872                     % Calculations for density, temperature modulations by magnetic ripple
0873                     %
0874                     cta = cos(theta)*calpha + sin(theta)*salpha;%for (Bx,By) -> (Brho,Bs) calculations
0875                     sta = -cos(theta)*salpha + sin(theta)*calpha;%for (Bx,By) -> (Brho,Bs) calculations
0876                     %
0877                     dctadrho = cos(theta)*dcalphadrho + sin(theta)*dsalphadrho;
0878                     dctadtheta = -sta + cos(theta)*dcalphadtheta + sin(theta)*dsalphadtheta;
0879                     %
0880                     dstadrho = -cos(theta)*dsalphadrho + sin(theta)*dcalphadrho;
0881                     dstadtheta = cta - cos(theta)*dsalphadtheta + sin(theta)*dcalphadtheta;
0882                     %
0883                     Brho0 = B*(cta*w_fit{1} + sta*w_fit{2});
0884                     dBrho0drho = dBdrho*(cta*w_fit{1} + sta*w_fit{2}) + B*(dctadrho*w_fit{1} + cta*dwdrho_fit{1} + dstadrho*w_fit{2} + sta*dwdrho_fit{2});
0885                     dBrho0dtheta = dBdtheta*(cta*w_fit{1} + sta*w_fit{2}) + B*(dctadtheta*w_fit{1} + cta*dwdtheta_fit{1} + dstadtheta*w_fit{2} + sta*dwdtheta_fit{2});
0886                     %
0887                     lambda = Brho0*gradrho*equil_fit.Rp/BPHI/nmin;
0888                     dlambdadrho = (dBrho0drho*gradrho + Brho0*dgradrhodrho - Brho0*gradrho*dBPHIdrho/BPHI)*equil_fit.Rp/BPHI/nmin;
0889                     dlambdadtheta = (dBrho0dtheta*gradrho + Brho0*dgradrhodtheta - Brho0*gradrho*dBPHIdtheta/BPHI)*equil_fit.Rp/BPHI/nmin;
0890                     %
0891                     rho_rip = rho + lambda*cos(nmin*phi); 
0892                     drhodrho_rip = 1 + dlambdadrho*cos(nmin*phi); 
0893                     drhodtheta_rip = dlambdadtheta*cos(nmin*phi);
0894                     drhodphi_rip = -nmin*lambda*sin(nmin*phi);
0895                     %
0896                     % Electron density
0897                     %
0898                     ne_sep = ppval_yp(equil_fit.ne_fit.pp_f,1);%value at the separatrix
0899                     dnedrho_sep = ppval_yp(equil_fit.ne_fit.pp_dfdrho,1);
0900                     %
0901                     if dnedrho_sep == 0,
0902                         rho_neeq0 = Inf;
0903                     else 
0904                         rho_neeq0 = - ne_sep/dnedrho_sep + 1;
0905                     end
0906                     %
0907                     if rho_rip <= 1,
0908                         ne_rip = ppval_yp(equil_fit.ne_fit.pp_f,rho_rip);%new value of the electron density that accounts for the flux surface modulation inside the separatrix
0909                         dnedrho_rip = ppval_yp(equil_fit.ne_fit.pp_dfdrho,rho_rip)*drhodrho_rip;
0910                         dnedtheta_rip = ppval_yp(equil_fit.ne_fit.pp_dfdrho,rho_rip)*drhodtheta_rip;
0911                         dnedphi_rip = ppval_yp(equil_fit.ne_fit.pp_dfdrho,rho_rip)*drhodphi_rip;
0912                     elseif (rho_rip > 1) & (rho_rip <= rho_neeq0),
0913                         ne_rip = ne_sep + dnedrho_sep*(rho_rip - 1);%SOL
0914                         dnedrho_rip = dnedrho_sep*drhodrho_rip;
0915                         dnedtheta_rip = dnedrho_sep*drhodtheta_rip;
0916                         dnedphi_rip = dnedrho_sep*drhodphi_rip;
0917                     else
0918                         ne_rip = 0;%vacuum
0919                         dnedrho_rip = 0;
0920                         dnedtheta_rip = 0;
0921                         dnedphi_rip = 0;
0922                     end
0923                     %
0924                     ne_tilde = ne_rip - ne;
0925                     dnedrho_tilde = dnedrho_rip - dnedrho;
0926                     dnedtheta_tilde = dnedtheta_rip;
0927                     dnedphi_tilde = dnedphi_rip;
0928                     %
0929                     ner_tilde = ne_tilde/ne;
0930                     dnerdrho_tilde = dnedrho_tilde/ne - ne_tilde*dnedrho/ne/ne;
0931                     dnerdtheta_tilde = dnedtheta_tilde/ne;
0932                     dnerdphi_tilde = dnedphi_tilde/ne;
0933                     %
0934                     fluctval.ner = fluctval.ner + ner_tilde;
0935                     fluctval.dnerdrho = fluctval.dnerdrho + dnerdrho_tilde;
0936                     fluctval.dnerdtheta = fluctval.dnerdtheta + dnerdtheta_tilde;
0937                     fluctval.dnerdphi = fluctval.dnerdphi + dnerdphi_tilde;
0938                     %
0939                     % Electron temperature
0940                     %
0941                     Te_sep = ppval_yp(equil_fit.Te_fit.pp_f,1);%value at the separatrix
0942                     dTedrho_sep = ppval_yp(equil_fit.Te_fit.pp_dfdrho,1);
0943                     %
0944                     if dTedrho_sep == 0,
0945                         rho_Teeq0 = Inf;
0946                     else 
0947                         rho_Teeq0 = - Te_sep/dTedrho_sep + 1;
0948                     end
0949                     %
0950                     if rho_rip <= 1,
0951                         Te_rip = ppval_yp(equil_fit.Te_fit.pp_f,rho_rip);%new value of the electron temperature that accounts for the flux surface modulation inside the separatrix
0952                         dTedrho_rip = ppval_yp(equil_fit.Te_fit.pp_dfdrho,rho_rip)*drhodrho_rip;
0953                         dTedtheta_rip = ppval_yp(equil_fit.Te_fit.pp_dfdrho,rho_rip)*drhodtheta_rip;
0954                         dTedphi_rip = ppval_yp(equil_fit.Te_fit.pp_dfdrho,rho_rip)*drhodphi_rip;
0955                     elseif (rho_rip > 1) & (rho_rip <= rho_Teeq0),
0956                         Te_rip = Te_sep + dTedrho_sep*(rho_rip - 1);%SOL
0957                         dTedrho_rip = dTedrho_sep*drhodrho_rip;
0958                         dTedtheta_rip = dTedrho_sep*drhodtheta_rip;
0959                         dTedphi_rip = dTedrho_sep*drhodphi_rip;
0960                     else
0961                         Te_rip = 0;%vacuum
0962                         dTedrho_rip = 0;
0963                         dTedtheta_rip = 0;
0964                         dTedphi_rip = 0;
0965                     end
0966                     %
0967                     Te_tilde = Te_rip - Te;
0968                     dTedrho_tilde = dTedrho_rip - dTedrho;
0969                     dTedtheta_tilde = dTedtheta_rip;
0970                     dTedphi_tilde = dTedphi_rip;
0971                     %
0972                     Ter_tilde = Te_tilde/Te;
0973                     dTerdrho_tilde = dTedrho_tilde/Te - Te_tilde*dTedrho/Te/Te;
0974                     dTerdtheta_tilde = dTedtheta_tilde/Te;
0975                     dTerdphi_tilde = dTedphi_tilde/Te;
0976                     %
0977                     fluctval.Ter = fluctval.Ter + Ter_tilde;
0978                     fluctval.dTerdrho = fluctval.dTerdrho + dTerdrho_tilde;
0979                     fluctval.dTerdtheta = fluctval.dTerdtheta + dTerdtheta_tilde;
0980                     fluctval.dTerdphi = fluctval.dTerdphi + dTerdphi_tilde;
0981                     %
0982                     % ion densities
0983                     %
0984                     zni_sep = ppval_yp(equil_fit.zni_fit.pp_f,1);%value at the separatrix
0985                     dznidrho_sep = ppval_yp(equil_fit.zni_fit.pp_dfdrho,1);
0986                     %
0987                     for ii = 1:length(equil_fit.zZi),
0988                         if dznidrho_sep(ii) == 0,
0989                             rho_znieq0(ii) = Inf;
0990                         else 
0991                             rho_znieq0(ii) = - zni_sep(ii)/dznidrho_sep(ii) + 1;
0992                         end
0993                     end
0994                     %
0995                     if rho_rip <= 1,
0996                         zni_rip = ppval_yp(equil_fit.zni_fit.pp_f,rho_rip);%new values of the ion densities that accounts for the flux surface modulation inside the separatrix
0997                         dznidrho_rip = ppval_yp(equil_fit.zni_fit.pp_dfdrho,rho_rip)*drhodrho_rip;
0998                         dznidtheta_rip = ppval_yp(equil_fit.zni_fit.pp_dfdrho,rho_rip)*drhodtheta_rip;
0999                         dznidphi_rip = ppval_yp(equil_fit.zni_fit.pp_dfdrho,rho_rip)*drhodphi_rip;
1000                     end
1001                     %
1002                     for ii = 1:length(equil_fit.zZi),    
1003                         if (rho_rip > 1) & (rho_rip <= rho_znieq0(ii)),
1004                             zni_rip(ii) = zni_sep(ii) + dznidrho_sep(ii)*(rho_rip - 1);%SOL
1005                             dznidrho_rip(ii) = dznidrho_sep(ii)*drhodrho_rip;
1006                             dznidtheta_rip(ii) = dznidrho_sep(ii)*drhodtheta_rip;
1007                             dznidphi_rip(ii) = dznidrho_sep(ii)*drhodphi_rip;
1008                         elseif (rho_rip > rho_znieq0(ii)),
1009                             zni_rip(ii) = 0;%vacuum
1010                             dznidrho_rip(ii) = 0;
1011                             dznidtheta_rip(ii) = 0;
1012                             dznidphi_rip(ii) = 0;
1013                         end
1014                     end
1015                     %
1016                     zni_tilde = zni_rip - zni;
1017                     dznidrho_tilde = dznidrho_rip - dznidrho;
1018                     dznidtheta_tilde = dznidtheta_rip;
1019                     dznidphi_tilde = dznidphi_rip;
1020                     %
1021                     znir_tilde = zni_tilde./zni;
1022                     dznirdrho_tilde = dznidrho_tilde./zni - zni_tilde.*dznidrho./zni./zni;
1023                     dznirdtheta_tilde = dznidtheta_tilde./zni;
1024                     dznirdphi_tilde = dznidphi_tilde./zni;
1025                     %
1026                     znir_tilde(isnan(znir_tilde)) = 0;
1027                     dznirdrho_tilde(isnan(dznirdrho_tilde)) = 0;
1028                     dznirdtheta_tilde(isnan(dznirdtheta_tilde)) = 0;
1029                     dznirdphi_tilde(isnan(dznirdphi_tilde)) = 0;
1030                     %
1031                     fluctval.znir = fluctval.znir + znir_tilde';
1032                     fluctval.dznirdrho = fluctval.dznirdrho + dznirdrho_tilde';
1033                     fluctval.dznirdtheta = fluctval.dznirdtheta + dznirdtheta_tilde';
1034                     fluctval.dznirdphi = fluctval.dznirdphi + dznirdphi_tilde';
1035                     %
1036                     % ion temperatures
1037                     %
1038                     zTi_sep = ppval_yp(equil_fit.zTi_fit.pp_f,1);%value at the separatrix
1039                     dzTidrho_sep = ppval_yp(equil_fit.zTi_fit.pp_dfdrho,1);
1040                     %
1041                     for ii = 1:length(equil_fit.zZi),
1042                         if dzTidrho_sep(ii) == 0,
1043                             rho_zTieq0(ii) = Inf;
1044                         else 
1045                             rho_zTieq0(ii) = - zTi_sep(ii)/dzTidrho_sep(ii) + 1;
1046                         end
1047                     end
1048                     %
1049                     if rho_rip <= 1,
1050                         zTi_rip = ppval_yp(equil_fit.zTi_fit.pp_f,rho_rip);%new values of the ion temperatures that accounts for the flux surface modulation inside the separatrix
1051                         dzTidrho_rip = ppval_yp(equil_fit.zTi_fit.pp_dfdrho,rho_rip)*drhodrho_rip;
1052                         dzTidtheta_rip = ppval_yp(equil_fit.zTi_fit.pp_dfdrho,rho_rip)*drhodtheta_rip;
1053                         dzTidphi_rip = ppval_yp(equil_fit.zTi_fit.pp_dfdrho,rho_rip)*drhodphi_rip;
1054                     end
1055                     %
1056                     for ii = 1:length(equil_fit.zZi),    
1057                         if (rho_rip > 1) & (rho_rip <= rho_zTieq0(ii)),
1058                             zTi_rip(ii) = zTi_sep(ii) + dzTidrho_sep(ii)*(rho_rip - 1);%SOL
1059                             dzTidrho_rip(ii) = dzTidrho_sep(ii)*drhodrho_rip;
1060                             dzTidtheta_rip(ii) = dzTidrho_sep(ii)*drhodtheta_rip;
1061                             dzTidphi_rip(ii) = dzTidrho_sep(ii)*drhodphi_rip;
1062                         elseif (rho_rip > rho_zTieq0(ii)),
1063                             zTi_rip(ii) = 0;%vacuum
1064                             dzTidrho_rip(ii) = 0;
1065                             dzTidtheta_rip(ii) = 0;
1066                             dzTidphi_rip(ii) = 0;
1067                         end
1068                     end
1069                     %
1070                     zTi_tilde = zTi_rip - zTi;
1071                     dzTidrho_tilde = dzTidrho_rip - dzTidrho;
1072                     dzTidtheta_tilde = dzTidtheta_rip;
1073                     dzTidphi_tilde = dzTidphi_rip;
1074                     %
1075                     zTir_tilde = zTi_tilde./zTi;
1076                     dzTirdrho_tilde = dzTidrho_tilde./zTi - zTi_tilde.*dzTidrho./zTi./zTi;
1077                     dzTirdtheta_tilde = dzTidtheta_tilde./zTi;
1078                     dzTirdphi_tilde = dzTidphi_tilde./zTi;
1079                     %
1080                     zTir_tilde(isnan(zTir_tilde)) = 0;
1081                     dzTirdrho_tilde(isnan(dzTirdrho_tilde)) = 0;
1082                     dzTirdtheta_tilde(isnan(dzTirdtheta_tilde)) = 0;
1083                     dzTirdphi_tilde(isnan(dzTirdphi_tilde)) = 0;
1084                     %
1085                     fluctval.zTir = fluctval.zTir + zTir_tilde';
1086                     fluctval.dzTirdrho = fluctval.dzTirdrho + dzTirdrho_tilde';
1087                     fluctval.dzTirdtheta = fluctval.dzTirdtheta + dzTirdtheta_tilde';
1088                     fluctval.dzTirdphi = fluctval.dzTirdphi + dzTirdphi_tilde';
1089                 end
1090             end
1091         end
1092     end
1093 end
1094 %
1095 
1096

Community support and wiki are available on Redmine. Last update: 18-Apr-2019.