other/Matlab/detection/plot_like_new.m

   1 %function plot_like_new
   2
   3         % function h=heat(T)
   4         c = 1;  % time constant (h) : how long the heat takes to drop to exp(-1)
   5         %********************************
   6         heat = @(T) (T<=0).*exp(min(T,0)/c);
   7         %********************************
   8         % T = fire arrival time
   9         % 0 if T > T_now = 0, 1 at T=0, exp decay for T>= 0
  10         figure(1);fplot(heat,[-5*c,c]);
  11         xlabel('time (h)');ylabel('heat fraction (1)')
  12
  13         % function p=prob(h)
  14         false_detection_rate = 0.001;  % = 1/(1+exp(b))
  15         b = log(1/false_detection_rate  -1);
  16         half_prob_heat = 0.05; % 1+exp(-a*h + b) = 2 for h=half_prob_heat
  17         a = b / half_prob_heat;
  18         %********************************
  19         prob = @(h) 1./(1+exp(-a*h + b));
  20         %********************************
  21         err1 = false_detection_rate - prob(0)
  22         err2 = 0.5 - prob(half_prob_heat)
  23         figure(2);fplot(prob,[0,half_prob_heat*5]);
  24         xlabel('heat fraction (1)');ylabel('probability of detection');
  25
  26         figure(3);fplot(@(T) prob(heat(-T)),[-c,10*c]);
  27         xlabel('time since fire arrival (h)');ylabel('probability of detection');
  28
  29         m = 10                  % mesh refinement
  30         dx=300/m, dy=300/m,    % grid step
  31         ix = 10*m, iy=10*m     % grid center indices
  32         nx = 2*ix-1, ny=2*iy-1     % grid dimension
  33         sigma = 500,     % geolocation error stdev
  34         R = 1            % rate of spread  m/s
  35         R = 3600*R       % rate of spread  m/h
  36
  37         [x,y]=ndgrid([0:nx-1]*dx,[0:ny-1]*dy);  % mesh
  38
  39         ssum = @(a) sum(a(:));    % utility: sum of array
  40
  41         % gaussian weights
  42         d2 = (x(ix,iy)-x).^2 + (y(ix,iy)-y).^2; % distance ^2
  43         w=exp(-d2/(sigma^2));                   % gaussian kernel
  44         w=w/sum(w(:));                          % weights add up to 1
  45
  46         %********************************
  47         tign = @(T)T+(x-x(ix,iy))/R;            % fire arrival time
  48         prbx = @(T)prob(heat(tign(T)));         % detection probability
  49         pgeo = @(T)ssum(w.*prob(heat(tign(T))));% weighted by geolocation error
  50         ageo = @(T)arrayfun(pgeo,T);            % same for array argument
  51         %********************************
  52         figure(4);fplot(@(T) log(ageo(-T)),[-c,10*c]);
  53         xlabel('time since fire arrival (h)');ylabel('probability of detection');
  54         %********************************
  55         %********************************
  56
  57
  58
  59
  60
  61  %end