1 function fuel_fraction=fuel_burnt(lfn,tign,tnow,fd,fuel_time)
\r
2 % lfn - level set function, size (2,2)
\r
3 % tign - time of ignition, size (2,2)
\r
5 % fd - mesh cell size (2)
\r
6 % fuel_time - time constant of fuel
\r
8 % output: approximation of
\r
11 % fuel_burnt = ---------- | ( 1 - exp( - ----------- ) ) dxdy
\r
12 % fd(1)*fd(2) | fuel_time
\r
17 % where R is the rectangle [0,fd(1)] * [0,fd(2)]
\r
18 % T is given by tign and L is given by lfn at the 4 vertices of R
\r
20 % note that fuel_left = 1 - fuel_burnt
\r
23 % exact in the case when T and L are linear
\r
24 % varies continuously with input
\r
25 % value of tign-tnow where lfn>0 ignored
\r
26 % assume T=0 when lfn=0
\r
28 for i=figs,figure(i),clf,hold off,end
\r
32 % nothing burning, nothing to do - most common case, put it first
\r
34 elseif all(lfn(:)<=0),
\r
36 % T=u(1)*x+u(2)*y+u(3)
\r
37 % set up least squares(A*u-v)*inv(C)*(A*u-v)->min
\r
42 v=tnow-tign(:) ; % time from ignition
\r
44 u = lscov(A,v,rw); % solve least squares to get coeffs of T
\r
45 residual=norm(A*u-v); % zero if T is linear
\r
47 display('fuel_burnt coefficients of u')
\r
49 uu = -u / fuel_time;
\r
50 fuel_fraction = 1 - exp(uu(3)) * intexp(uu(1)*fd(1)) * intexp(uu(2)*fd(2));
\r
51 if(fuel_fraction<0 | fuel_fraction>1),
\r
52 warning('fuel_fraction should be between 0 and 1')
\r
56 % part of cell is burning - the interesting case
\r
57 % set up a list of points with the corresponding values of T and L
\r
58 xylist=zeros(8,2); % allocate list of points
\r
61 xx=[-fd(1) fd(1) fd(1) -fd(1) -fd(1)]/2; % make center (0,0)
\r
62 yy=[-fd(2) -fd(2) fd(2) fd(2) -fd(2)]/2; % cyclic, counterclockwise
\r
63 ii=[1 2 2 1 1]; % indices of corners, cyclic, counterclockwise
\r
67 lfn0=lfn(ii(k),jj(k));
\r
68 lfn1=lfn(ii(k+1),jj(k+1));
\r
71 xylist(points,:)=[xx(k),yy(k)]; % add corner to list
\r
72 Tlist(points)=tnow-tign(ii(k),jj(k)); % time since ignition
\r
77 % coordinates of intersection of fire line with segment k k+1
\r
78 %lfn(t)=lfn0 + t*(lfn1-lfn0)=0
\r
80 x0=xx(k)+(xx(k+1)-xx(k))*t;
\r
81 y0=yy(k)+(yy(k+1)-yy(k))*t;
\r
82 xylist(points,:)=[x0,y0];
\r
83 Tlist(points)=0; % now at ignition
\r
84 Llist(points)=0; % at fireline
\r
87 % make the lists circular and trim to size
\r
88 Tlist(points+1)=Tlist(1);
\r
89 Tlist=Tlist(1:points+1);
\r
90 Llist(points+1)=Llist(1);
\r
91 Llist=Llist(1:points+1);
\r
92 xylist(points+1,:)=xylist(1,:);
\r
93 xylist=xylist(1:points+1,:);
\r
94 for k=1:5,lfnk(k)=lfn(ii(k),jj(k));end
\r
96 plot3(xx,yy,lfnk,'k')
\r
98 plot3(xylist(:,1),xylist(:,2),Tlist,'m--o')
\r
99 plot3(xylist(:,1),xylist(:,2),Llist,'g--o')
\r
101 plot(xylist(:,1),xylist(:,2),'-.r+')
\r
102 legend('lfn on whole cell','tign-tnow on burned area',...
\r
103 'lfn on burned area', 'cell boundary','burned area boundary')
\r
104 patch(xylist(:,1),xylist(:,2),zeros(points+1,1),250,'FaceAlpha',0.3)
\r
105 patch(xylist(:,1),xylist(:,2),Tlist,100,'FaceAlpha',0.3)
\r
106 hold off,grid on,drawnow,pause(0.5)
\r
107 % set up least squares
\r
108 A=[xylist(1:points,1:2),ones(points,1)];
\r
113 dist(j)=norm(xylist(i,:)-xylist(j,:));
\r
115 dist(j)=max(fd); % large
\r
118 rw(i)=1+min(dist); % weight = 1+min dist from other nodes
\r
120 u = lscov(A,v,rw); % solve least squares to get coeffs of T
\r
122 display('fuel_burnt coefficients of u')
\r
123 residual=norm(A*u-v); % should be zero if T and L are linear
\r
124 nr=sqrt(u(1)*u(1)+u(2)*u(2));
\r
127 Q=[c, s; -s, c]; % rotation such that Q*u(1:2)=[something;0]
\r
129 errQ=ut(2); % should be zero
\r
130 ae=-ut(1)/fuel_time;
\r
131 ce=-u(3)/fuel_time; % -T(xt,yt)/fuel_time=ae*xt+ce
\r
132 cet=ce; %keep ce for later
\r
133 xytlist=xylist*Q'; % rotate the points in the list
\r
134 xt=sort(xytlist(1:points,1)); % sort ascending in x
\r
137 % integrate the vertical slice from xt1 to xt2
\r
139 plot(xytlist(:,1),xytlist(:,2),'-o'),grid on,hold on
\r
142 plot([xt1,xt1],sum(fd)*[-1,1]/3,'y')
\r
143 plot([xt2,xt2],sum(fd)*[-1,1]/3,'y')
\r
144 if(xt2-xt1>100*eps*max(fd)), % slice of nonzero width
\r
145 % find slice height as h=ah*x+ch
\r
148 for s=1:points % pass counterclockwise
\r
149 xts=xytlist(s,1); % start point of the line
\r
151 xte=xytlist(s+1,1); % end point of the line
\r
152 yte=xytlist(s+1,2);
\r
153 if (xts>xt1 & xte > xt1) | (xts<xt2 & xte < xt2)
\r
154 plot([xts,xte],[yts,yte],'--k')
\r
155 % that is not the one
\r
156 else % line y=a*x+c through (xts,yts), (xte,yte)
\r
157 a=(yts-yte)/(xts-xte);
\r
158 c=(xts*yte-xte*yts)/(xts-xte);
\r
159 if xte<xts % upper boundary
\r
162 plot([xts,xte],[yts,yte],'g')
\r
166 else % lower boundary
\r
169 plot([xts,xte],[yts,yte],'m')
\r
174 if(lower~=1|upper~=1),
\r
175 error('slice does not have one upper and one lower line')
\r
177 ce=cet;% use kept ce
\r
178 % shift samll amounts to avoid negative fuel burnt
\r
179 if ae*xt1+ce > 0;%disp('ae*xt1+ce');disp(ae*xt1+ce);pause;
\r
180 ce=ce-(ae*xt1+ce);end;
\r
181 if ae*xt2+ce > 0;%disp('ae*xt2+ce');disp(ae*xt2+ce);pause;
\r
182 ce=ce-(ae*xt2+ce);end;
\r
187 patch([xt1,xt2,xt2,xt1],...
\r
188 [alow*[xt1,xt2]+clow,aupp*[xt2,xt1]+cupp,],k*10)
\r
189 axis equal,drawnow, pause(0.5)
\r
191 x=[xt1:(xt2-xt1)/100:xt2];
\r
192 plot(x,1-exp(ae*x+ce),x,ah*x+ch,x,...
\r
193 (ah*x+ch).*(1-exp(ae*x+ce)),[xt1,xt2],[0,0],'+k')
\r
194 xlabel x,legend('burned frac','slice height','height*burned','dividers')
\r
195 hold on,grid on,drawnow,pause(0.5)
\r
196 % integrate (ah*x+ch)*(1-exp(ae*x+ce) from xt1 to xt2
\r
197 % numerically sound for ae->0, ae -> infty
\r
198 % this can be important for different model scales
\r
199 % esp. if someone runs the model in single precision!!
\r
200 % s1=int((ah*x+ch),x,xt1,xt2)
\r
201 s1=(xt2-xt1)*((1/2)*ah*(xt2+xt1)+ch);
\r
202 % s2=int((ch)*(-exp(ae*x+ce)),x,xt1,xt2)
\r
204 s2=-ch*exp(ae*(xt1+ceae))*(xt2-xt1)*intexp(ae*(xt2-xt1));
\r
205 % s3=int((ah*x)*(-exp(ae*x+ce)),x,xt1,xt2)
\r
206 % s3=int((ah*x)*(-exp(ae*(x+ceae))),x,xt1,xt2)
\r
207 % expand in Taylor series around ae=0
\r
208 % collect(expand(taylor(int(x*(-exp(ae*(x+ceae))),x,xt1,xt2)*ae^2,ae,4)/ae^2),ae)
\r
209 % =(1/8*xt1^4+1/3*xt1^3*ceae+1/4*xt1^2*ceae^2-1/8*xt2^4-1/3*xt2^3*ceae-1/4*xt2^2*ceae^2)*ae^2
\r
210 % + (-1/3*xt2^3-1/2*xt2^2*ceae+1/3*xt1^3+1/2*xt1^2*ceae)*ae
\r
211 % + 1/2*xt1^2-1/2*xt2^2
\r
213 % coefficient at ae^2 in the expansion, after some algebra
\r
214 a2=(xt1-xt2)*((1/4)*(xt1+xt2)*ceae^2+(1/3)*(xt1^2+xt1*xt2+xt2^2)*ceae+(1/8)*(xt1^3+xt1*xt2^2+xt1^2*xt2+xt2^3));
\r
217 % since ae*xt1+ce<=0 ae*xt2+ce<=0 all fine for large ae
\r
218 % for ae, ce -> 0 rounding error approx eps/ae^2
\r
219 s3=(exp(ae*(xt1+ceae))*(ae*xt1-1)-exp(ae*(xt2+ceae))*(ae*xt2-1))/ae^2;
\r
220 % we do not worry about rounding as xt1 -> xt2, then s3 -> 0
\r
222 % coefficient at ae^1 in the expansion
\r
223 a1=(xt1-xt2)*((1/2)*ceae*(xt1+xt2)+(1/3)*(xt1^2+xt1*xt2+xt2^2));
\r
224 % coefficient at ae^0 in the expansion for ae->0
\r
225 a0=(1/2)*(xt1-xt2)*(xt1+xt2);
\r
226 s3=a0+a1*ae+a2*ae^2; % approximate the integral
\r
229 fuel_fraction=fuel_fraction+s1+s2+s3;
\r
230 if(fuel_fraction<0 | fuel_fraction>(fd(1)*fd(2))),
\r
231 fuel_fraction,(fd(1)*fd(2)),s1,s2,s3
\r
232 warning('fuel_fraction should be between 0 and 1')
\r
237 fuel_fraction=fuel_fraction/(fd(1)*fd(2));
\r
238 display('fuel_burnt2')
\r
239 %fuel_fraction=1-fuel_fraction;
\r
241 for i=figs,figure(i),hold off,end
\r
244 function s=intexp(ab)
\r
245 % function s=intexp(ab)
\r
246 % s = (1/b)*int(exp(a*x),x,0,b) ab = a*b
\r
248 if eps < abs(ab)^3/6,
\r
249 s = (exp(ab)-1)/(ab); % rounding error approx eps/(a*b) for small a*b
\r
251 s = 1 + ab/2+ab^2/6; % last term (a*b)^2/6 for small a*b
\r