cycling/multi_make.m

   1 function mt = multi_make(ps,alpha,p_param,grid_fraction)
   2
   3 %do make_tign on set of grids of decreasing size
   4 gs = [2000,1000,750,600,500];
   5 %compute grid fractions gf = gs/dlon
   6 red = ps.red;
   7 E = wgs84Ellipsoid;
   8 dlon= distance(red.min_lat,red.min_lon,red.min_lat,red.max_lon,E);
   9 dlat= distance(red.min_lat,red.min_lon,red.max_lat,red.min_lon,E);
  10 [m,n] = size(red.tign);
  11 gf = (dlon/m)/gs(1);
  12
  13 %make series of estimates
  14 grids = length(gs);
  15 %10 steps, max
  16 grids = 10;
  17 re1 = 1;
  18 do_plots = 0;
  19 for i = 1:grids
  20     [t,e] = make_tign(ps,alpha,p_param,gf);
  21     if do_plots == 1
  22        close all
  23        tstr = sprintf('Iteration %d',i);
  24        savestr = sprintf('multi_%d',i);
  25        %limit plot region
  26        t_max = max(t(:));
  27        msk = t<t_max;
  28        min_lon = min(ps.red.fxlong(msk));
  29        max_lon = max(ps.red.fxlong(msk));
  30        lon_buff = 0.2*(max_lon-min_lon);
  31        min_lat = min(ps.red.fxlat(msk));
  32        max_lat = max(ps.red.fxlat(msk));
  33        lat_buff = 0.2*(max_lat-min_lat);
  34        figure(312)
  35        mesh(ps.red.fxlong(1:5:end,1:5:end),ps.red.fxlat(1:5:end,1:5:end),t(1:5:end,1:5:end));
  36        hold on
  37        scatter3(ps.points(:,2),ps.points(:,1),ps.points(:,3),'r*');
  38        xlim([min_lon-lon_buff max_lon+lon_buff]);
  39        ylim([min_lat-lat_buff max_lat+lat_buff]);
  40        xlabel('Lon (deg)');
  41        ylabel('Lat (deg)');
  42        zlabel('Time (datenum)');
  43        title(tstr);
  44        savefig(savestr);
  45        saveas(gcf,[savestr '.png']);
  46        hold off
  47     end
  48     ps.red.tign = t;
  49     mt(i).t = t;
  50     mt(i).e = e;
  51     if i < grids
  52         %gf = gf*gs(i)/gs(i+1);
  53         gf = 4/3*gf;
  54         close all
  55     end
  56     if i > 1
  57         re1 = norm(mt(i).t-mt(i-1).t)/norm(mt(i-1).t-min(mt(i-1).t(:)));
  58     end
  59     fprintf('relative error from previous estimate : \n',re1)
  60     if re1 < 0.001
  61         fprintf('not changing.....  exit loop  \n')
  62         break
  63     end
  64     if i == grids
  65         fprintf('Max steps performed\n')
  66     end
  67 end
  68
  69
  70 %fix the ripple on the top
  71 r_mask = mt(end).t >=max(mt(end).t(:))-0.05;
  72 mt(end).t(r_mask) = max(mt(end).t(:));
  73
  74 %compute relative error
  75 re = norm(red.tign-mt(end).t)/(norm(red.tign-red.start_datenum));
  76 fprintf('Relative error: %f\n',re)
  77
  78 end % function