delta_wave_optim_jd

0001 function [delta,ntargetsh] = delta_wave_optim_jd(wave,equil_fit,ntargets)
0002 %
0003 ray = wave.rays{1};
0004 if ray.P0_2piRp == ray.sP_2piRp_lin(end),% no absorption
0005     ray.sP_2piRp_lin(end) = ray.P0_2piRp*(1 - 10*eps);
0006 end
0007 Ph_2piRp_lin = (ray.P0_2piRp + ray.sP_2piRp_lin(end))/2;
0008 %
0009 smask = 1 + [0,find(diff(ray.sP_2piRp_lin) < 0)];% monotonic evolution
0010 %
0011 rhoh = interp1(ray.sP_2piRp_lin(smask),sqrt(ray.spsin(smask)),Ph_2piRp_lin);
0012 thetah = interp1(ray.sP_2piRp_lin(smask),ray.stheta(smask),Ph_2piRp_lin);
0013 nparh = interp1(ray.sP_2piRp_lin(smask),ray.sNpar(smask),Ph_2piRp_lin);
0014 %
0015 % xh = interp1(ray.sP_2piRp_lin(smask),ray.sx(smask),Ph_2piRp_lin);
0016 % yh = interp1(ray.sP_2piRp_lin(smask),ray.sy(smask),Ph_2piRp_lin);
0017 %
0018 x_a0h = ppval_yp(equil_fit.x_fit.pp_a0,rhoh);
0019 x_anh = ppval_yp(equil_fit.x_fit.pp_an,rhoh);
0020 x_bnh = ppval_yp(equil_fit.x_fit.pp_bn,rhoh);
0021 %
0022 y_a0h = ppval_yp(equil_fit.y_fit.pp_a0,rhoh);
0023 y_anh = ppval_yp(equil_fit.y_fit.pp_an,rhoh);
0024 y_bnh = ppval_yp(equil_fit.y_fit.pp_bn,rhoh);
0025 %
0026 xh = x_a0h + x_anh*cos(equil_fit.x_fit.n*thetah) + x_bnh*sin(equil_fit.x_fit.n*thetah);
0027 yh = y_a0h + y_anh*cos(equil_fit.y_fit.n*thetah) + y_bnh*sin(equil_fit.y_fit.n*thetah);
0028 %
0029 ntargetsh(1) = rhoh;
0030 ntargetsh(2) = thetah;
0031 ntargetsh(3) = nparh;
0032 %
0033 if ~isempty(ntargets),
0034     %
0035     if ~isnan(ntargets(1)),
0036         rho = ntargets(1);
0037     else
0038         rho = rhoh;
0039     end
0040     %
0041     if ~isnan(ntargets(2)),
0042         theta = ntargets(2);
0043     else
0044         theta = thetah;
0045     end
0046     %
0047     if ~isnan(ntargets(3)),
0048         npar = ntargets(3);
0049     else
0050         npar = nparh;
0051     end
0052     %
0053     x_a0 = ppval_yp(equil_fit.x_fit.pp_a0,rho);
0054     x_an = ppval_yp(equil_fit.x_fit.pp_an,rho);
0055     x_bn = ppval_yp(equil_fit.x_fit.pp_bn,rho);
0056     %
0057     y_a0 = ppval_yp(equil_fit.y_fit.pp_a0,rho);
0058     y_an = ppval_yp(equil_fit.y_fit.pp_an,rho);
0059     y_bn = ppval_yp(equil_fit.y_fit.pp_bn,rho);
0060     %
0061     x = x_a0 + x_an*cos(equil_fit.x_fit.n*theta) + x_bn*sin(equil_fit.x_fit.n*theta);
0062     y = y_a0 + y_an*cos(equil_fit.y_fit.n*theta) + y_bn*sin(equil_fit.y_fit.n*theta);
0063     %
0064     delta = sqrt((x - xh).^2 + (y - yh).^2 + (npar - nparh).^2);
0065     %
0066 else
0067     delta = [];
0068 end