makegrid_dke_jd

 This function generates a non-uniform grid.

   INPUTS:

       - xmin: initial grid position
       - nmin: minimum number of fine grid points near xmin
       - dxmin: minimum grid step near xmin
       - xmax: final grid position 
       - nmax: minimum number of fine grid points near xmax
       - dxmax: minimum grid step near xmax
       - ntot: minimum total number of points, including ends 
       - dtop: maximum grid step
       - grid_power: maximum ratio between size of two consecutive steps
       (y_power > sqrt(2) encouraged)

   OUTPUTS:

       - x: non-uniform grid

  by Joan Decker <jodecker@alum.mit.edu> (MIT/RLE) and Yves Peysson <yves.peysson@cea.fr> (CEA/DRFC)

0001 function x = makegrid_dke_jd(xmin,nmin,dxmin,xmax,nmax,dxmax,ntot,dtop,grid_power);
0002 %
0003 % This function generates a non-uniform grid.
0004 %
0005 %   INPUTS:
0006 %
0007 %       - xmin: initial grid position
0008 %       - nmin: minimum number of fine grid points near xmin
0009 %       - dxmin: minimum grid step near xmin
0010 %       - xmax: final grid position
0011 %       - nmax: minimum number of fine grid points near xmax
0012 %       - dxmax: minimum grid step near xmax
0013 %       - ntot: minimum total number of points, including ends
0014 %       - dtop: maximum grid step
0015 %       - grid_power: maximum ratio between size of two consecutive steps
0016 %       (y_power > sqrt(2) encouraged)
0017 %
0018 %   OUTPUTS:
0019 %
0020 %       - x: non-uniform grid
0021 %
0022 %  by Joan Decker <jodecker@alum.mit.edu> (MIT/RLE) and Yves Peysson <yves.peysson@cea.fr> (CEA/DRFC)
0023 %
0024 if nargin < 9
0025     error('not enough arguments')
0026 end
0027 %
0028 if ntot > nmin + nmax
0029     warning('ntot = nmin + nmax enforced')
0030     ntot = nmin + nmax;
0031 end
0032 %
0033 if abs((dxmax - dxmin)/dxmax) < 1e-10
0034     dxmax = dxmin;
0035 end
0036 %
0037 xtot = xmax - xmin;
0038 %
0039 if dxmax < dxmin
0040     nc = nmin;
0041     nmin = nmax;
0042     nmax = nc;
0043     %
0044     dxc = dxmin;
0045     dxmin = dxmax;
0046     dxmax = dxc;
0047     %
0048     clear nc dxc
0049     inv = 1;
0050 else
0051     inv = 0;
0052 end
0053 %
0054 grid_ok = 0;
0055 %
0056 xin = powerdiff_dke_jd(nmin*dxmin, dxmin, xtot - nmax*dxmax, dxmax, dtop, grid_power);
0057 if isempty(xin),
0058     while grid_ok == 0 & dxmax > dxmin  
0059         for imax = nmax:-1:0,
0060             xin = powerdiff_dke_jd(nmin*dxmin, dxmin, xtot - imax*dxmax, dxmax, dxmax, grid_power);
0061             if isempty(xin) | length(xin) < (ntot - nmin - imax + 1) | diff(xin(length(xin)-1:length(xin))) < dxmin, 
0062                 for imin = 1:(ntot - imax - nmin - 1),
0063                     xin = powerdiff_dke_jd((nmin + imin)*dxmin, dxmin, xtot - imax*dxmax, dxmax, dxmax, grid_power);
0064                     if length(xin) >= (ntot - imax - nmin - imin + 1) & diff(xin(length(xin)-1:length(xin))) >= dxmin,
0065                         nmin = nmin + imin;
0066                         grid_ok = 1;
0067                         break
0068                     end
0069                 end
0070             else 
0071                 grid_ok = 1;
0072             end
0073             if grid_ok == 1,
0074                 nmax = imax;
0075                 break
0076             end
0077         end
0078         if grid_ok == 0,
0079             dxmax = max([dxmin,dxmax/grid_power]);
0080         end
0081     end
0082     if grid_ok == 0,
0083         nmin = 1;
0084         nmax = 1;
0085         xintot = xtot - nmin*dxmin - nmax*dxmin;
0086         itot = max([ntot - nmin - nmax,ceil(xintot/dxmin),1]);
0087         dxin = xintot/itot;
0088         if dxin > dxmin/grid_power
0089             xin = linspace(nmin*dxmin,xtot - nmax*dxmin,itot + 1);
0090         else
0091             nmin = 0;
0092             nmax = 0;
0093             ntot = max([ntot,ceil(xtot/dxmin)]);
0094             xin = linspace(0,xtot,ntot + 1);
0095         end
0096     end
0097 end
0098 %
0099 rampmin = (0: 1:nmin-1)*dxmin;
0100 rampmax = xtot - (nmax-1:-1:0)*dxmax;
0101 %
0102 x = [rampmin,xin,rampmax];  
0103 dx = diff(x);
0104 %
0105 if inv == 0;
0106     x = cumsum([xmin,dx]);
0107 else
0108     x = cumsum([xmin,fliplr(dx)]);
0109 end
0110 %
0111 x(length(x)) = xmax;%to compensate numerical errors from cumsum.