varspectrum_jd

 This function discretizes the spectrum characterized by the base
 Gaussian spectrum (npar0, dnpar, P0) and a Gaussian npar distribution
 characterized by dnpar0

 by J. Decker and Y. Peysson

0001 function [nP0n,nnpar0,bdnpar] = varspectrum_jd(P0,npar0,dnpar,dnpar0,snpar,fnpar,opt)
0002 %
0003 % This function discretizes the spectrum characterized by the base
0004 % Gaussian spectrum (npar0, dnpar, P0) and a Gaussian npar distribution
0005 % characterized by dnpar0
0006 %
0007 % by J. Decker and Y. Peysson
0008 %
0009 %
0010 if nargin < 7,
0011     opt = 1;
0012 end
0013 %
0014 if nargin < 6,
0015     fnpar = sqrt(log(100));
0016 end
0017 %
0018 if nargin < 5,
0019     snpar = 2*log(2);
0020 end
0021 %
0022 if nargin < 4,
0023     dnpar0 = 0.0;
0024 end
0025 %
0026 dnparb = sqrt(dnpar^2 + dnpar0^2);
0027 %
0028 if real(snpar) > 0,% snpar*dnpar is the separation between rays
0029     %
0030     %nparmin = npar0 - fnpar*dnparb;
0031     nparmax = npar0 + fnpar*dnparb;
0032     %
0033     %nmin = fix((npar0 - nparmin)/snpar/dnpar);
0034     nmax = fix((nparmax - npar0)/snpar/dnpar);
0035     %
0036 else % imag(snpar) is the number of rays
0037     %
0038     nparmax = npar0 + fnpar*dnpar0;
0039     %
0040     nmax = imag(snpar);
0041     snpar = (nparmax - npar0)/(nmax + 1/2)/dnpar;
0042     %
0043 end
0044 %
0045 n = -nmax:nmax;%-nmin:nmax;
0046 nn = length(n);
0047 %
0048 nnpar0 = npar0 + snpar*n*dnpar;
0049 %
0050 if opt == 0,
0051     %
0052     nnpar_S = [-Inf,(nnpar0(1:nn-1) + nnpar0(2:nn))/2,Inf];
0053     nnparb_S = (nnpar_S - npar0)/dnparb;
0054     nP0n = (erf(nnparb_S(2:end)) - erf(nnparb_S(1:end-1)))/2;
0055     %
0056 elseif opt == 1,
0057     %
0058     Xhnpar0 = repmat(nnpar0,[nn,1]);
0059     Xvnpar0 = repmat(nnpar0.',[1,nn]);
0060     XPn = exp(-(Xhnpar0 - Xvnpar0).^2/dnpar^2)/sqrt(pi)/dnpar;
0061     nPb = exp(-(nnpar0.' - npar0).^2/dnparb^2)/sqrt(pi)/dnparb;
0062     nP0n = (XPn\nPb).';
0063     %
0064     SP0n = sum(nP0n);
0065     disp(['- sum on P0n from inversion : ',num2str(SP0n)]);
0066     nP0n = nP0n/SP0n;
0067     %
0068 elseif opt == 2,
0069     %
0070     nP0n = exp(-(nnpar0 - npar0).^2/dnpar0^2)/sqrt(pi)/dnpar0;
0071     %
0072     SP0n = sum(nP0n);
0073     nP0n = nP0n/SP0n;
0074     %
0075 end
0076 %
0077 nP0n = nP0n*P0;
0078 bdnpar = repmat(dnpar,[1,nn]);
0079 %