radialdampingprofile_jd

 calculates the radial power deposition profile of a ray for which the
 damping is already calculated

 by J. Decker (joan.decker@cea.fr) and Y. Peysson (yves.peysson@cea.fr)

0001 function [sdPdrho_dep,rho_dep] = radialdampingprofile_jd(srho,sP,rho_dep_S)
0002 %
0003 % calculates the radial power deposition profile of a ray for which the
0004 % damping is already calculated
0005 %
0006 % by J. Decker (joan.decker@cea.fr) and Y. Peysson (yves.peysson@cea.fr)
0007 %
0008 if nargin < 3,
0009     rho_dep_S = 40;
0010 end
0011 %
0012 if nargin < 2,
0013     error('Not enough input arguments');
0014 end
0015 %
0016 if length(rho_dep_S) == 1,% number of points is given
0017     rho_dep_S = linspace(0,1,rho_dep_S+1);
0018 end
0019 %
0020 if rho_dep_S(1) > 0.000001,
0021     rho_dep_S = [0,rho_dep_S];
0022 end
0023 %
0024 if rho_dep_S(end) < 0.999999,
0025     rho_dep_S = [rho_dep_S,1];
0026 end
0027 %
0028 ns = length(srho);
0029 %
0030 rho_dep = (rho_dep_S(1:end-1) + rho_dep_S(2:end))/2;
0031 drho_dep = diff(rho_dep_S);
0032 nr_dep = length(rho_dep);
0033 %
0034 sdPdrho_dep = zeros(1,nr_dep);
0035 srhop = srho(1);
0036 irho_depp = find(srhop > rho_dep_S(1:end-1) & srhop <= rho_dep_S(2:end));
0037 for is = 2:ns,
0038     srhom = srhop;
0039     irho_depm = irho_depp;
0040     srhop = srho(is);
0041     irho_depp = find(srhop > rho_dep_S(1:end-1) & srhop <= rho_dep_S(2:end));
0042     %
0043     sPm = sP(is-1);
0044     sPp = sP(is);
0045     %
0046     if isempty(irho_depm) || isempty(irho_depp),
0047         continue
0048     elseif irho_depm == irho_depp,
0049         sdPdrho_dep(irho_depm) = sdPdrho_dep(irho_depm) + ...
0050             (sPm - sPp)/drho_dep(irho_depm);
0051     elseif srhop > srhom,
0052         sdPdrho_dep(irho_depm) = sdPdrho_dep(irho_depm) + ...
0053             (sPm - sPp)*(rho_dep_S(irho_depm + 1) - srhom)/(srhop - srhom)/drho_dep(irho_depm);
0054         sdPdrho_dep(irho_depm+1:irho_depp-1) = sdPdrho_dep(irho_depm+1:irho_depp-1) + ...
0055             (sPm - sPp)*(rho_dep_S(irho_depm+2:irho_depp) - rho_dep_S(irho_depm+1:irho_depp-1))/(srhop - srhom)./drho_dep(irho_depm+1:irho_depp-1);
0056         sdPdrho_dep(irho_depp) = sdPdrho_dep(irho_depp) + ...
0057             (sPm - sPp)*(srhop - rho_dep_S(irho_depp))/(srhop - srhom)/drho_dep(irho_depp);
0058     else
0059         sdPdrho_dep(irho_depm) = sdPdrho_dep(irho_depm) + ...
0060             (sPm - sPp)*(srhom - rho_dep_S(irho_depm))/(srhom - srhop)/drho_dep(irho_depm);
0061         sdPdrho_dep(irho_depp+1:irho_depm-1) = sdPdrho_dep(irho_depp+1:irho_depm-1) + ...
0062             (sPm - sPp)*(rho_dep_S(irho_depp+2:irho_depm) - rho_dep_S(irho_depp+1:irho_depm-1))/(srhom - srhop)./drho_dep(irho_depp+1:irho_depm-1);
0063         sdPdrho_dep(irho_depp) = sdPdrho_dep(irho_depp) + ...
0064             (sPm - sPp)*(rho_dep_S(irho_depp + 1) - srhop)/(srhom - srhop)/drho_dep(irho_depp);
0065     end     
0066 end
0067 %
0068