fuel_left/fuel_quad.m

   1 function fq = fuel_quad(samp,m,n)
   2 %computes quadrature
   3 %samp - matrix with [tign00 tign01 tign10 tign11 time_now fuel_time_cell]
   4 %or sp -- , mtrix with [tign00 tign01 tign10 tign11 fuel_time_cell]where
   5 %   t00 is processed so that t00 = time_now - t00
   6 %m - size of quadrature grid (m x m)
   7 %n - number of integrals to evaluate
   8
   9 proc = 1; %
  10 %proc = input_num('Is data processed? 1 = yes',0);
  11
  12 %transpose if samplles are from column vectors
  13 sz = size(samp);
  14 if sz(1)<sz(2)
  15     samp = samp';
  16 end
  17
  18
  19 %corners of cell [a,b]X[a,b]
  20 a = 0.0;
  21 b = 1.0;
  22 %grid spacing
  23 dx = (b-a)/m;
  24 dy = dx;
  25
  26 if ~exist('n','var')
  27     n = length(samp);
  28 end
  29 fq = zeros(1,n);
  30 %loops
  31 tic
  32 for o = 1:n % outer loop
  33     %fire parameters
  34     tign00 = samp(o,1);
  35     tign01 = samp(o,2);
  36     tign10 = samp(o,3);
  37     tign11 = samp(o,4);
  38     if proc == 1
  39         fuel_time_cell = samp(o,5);
  40     else
  41         time_now  = samp(o,5);
  42         fuel_time_cell = samp(o,6);
  43     end
  44     %integral sum --> s
  45     s = 0;
  46     %midpoint method
  47     for i = 1:m
  48         for j = 1:m
  49             x = a + dx/2.0 + (i-1.0)*dx;
  50             y = a + dy/2.0 + (j-1.0)*dy;
  51             % interpolate tign
  52             intt = (1.0-x)*( (1.0-y)*tign00 + y*tign01 ) + x*( (1.0-y)*tign10 + y*tign11 );
  53             % t is the time since fire arrived = time_now - intt
  54             if proc == 1
  55                 t = intt;
  56             else
  57                 t = time_now - intt;
  58             end
  59             if t >= 0
  60                 bf = (1 - exp(-t/fuel_time_cell));
  61                 % make sum update
  62                 s = s + bf;
  63             end %if
  64         end % for j
  65     end % for i
  66     %multiply differencials
  67     fq(1,o) = s*dx*dy;
  68
  69     %matlab integral2 method
  70     %intt = (1.0-x)*( (1.0-y)*tign00 + y*tign01 ) + x*( (1.0-y)*tign10 + y*tign11 );
  71
  72
  73 end % for o
  74 t = toc;
  75 %fprintf('Elepsed time is %f , Each quad was %f \n', toc, toc/n);