perimeter/perimeter.m

   1 function result=perimeter(long,lat,uf,vf,dzdxf,dzdyf,time_now,bound)
   2
   3 % Volodymyr Kondratenko           December 8 2012
   4
   5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   6
   7 % The function creates the initial matrix of times of ignitions
   8 %                       given the perimeter and time of ignition at the points of the perimeter and
   9 %                       is using wind and terrain gradient in its calculations
  10 % Input: We get it after reading the data using function read_file_perimeter.m
  11 %
  12 %        long                   FXLONG*UNIT_FXLONG, longtitude coordinates of the mesh converted into meters
  13 %        lat                    FXLAT*UNIT_FXLAT, latitude coordinates of the mesh converted into meters
  14 %        uf,vf                  horizontal wind velocity vectors of the points of the mesh
  15 %        dzdxf,dzdyf    terrain gradient of the points of the mesh
  16 %        time_now               time of ignition on the fireline (fire perimeter)
  17 %        bound                  set of ordered points of the fire perimeter 1st=last
  18 %                                               bound(i,1)-horisontal; bound(i,1)-vertical coordinate
  19 %
  20 % Output:   Matrix of time of ignition
  21 %
  22 % Code:
  23
  24 %addpath('../../other/Matlab/util1_jan');
  25 %addpath('../../other/Matlab/netcdf');
  26 tic
  27 fuels % This function is needed to create fuel variable, that contains all the characteristics
  28       % types of fuel, this function lies in the same folder where you run the main_function.m
  29           % (where is it originally located/can be created?)
  30
  31 format long
  32
  33 bnd_size=size(bound);
  34 n=size(long,1);
  35 m=size(long,2);
  36
  37 %tign=zeros(n,m);      % "time of ignition matrix" of the nodes
  38 %A=zeros(n,m);         % flag matrix of the nodes, A(i,j)=1 if the time of ignition of the
  39                       % point (i,j) was updated at least once
  40
  41 'started'
  42 %  IN - matrix, that shows, whether the point is inside (IN(x,y)>0) the burning region
  43 %  or outside (IN(x,y)<0)
  44 %  ON - matrix that, shows whether the point is on the boundary or not
  45 %  Both matrices evaluated using "polygon", coefficients are multiplied by
  46 %  10^6, because the function looses acuracy when it deals with decimals
  47
  48 xv=bound(:,1);
  49 yv=bound(:,2);
  50 xv=xv*100000;
  51 yv=yv*100000;
  52 lat1=lat*100000;
  53 long1=long*100000;
  54 [IN,ON] = inpolygon(long1,lat1,xv,yv);
  55
  56 % Code
  57
  58 [ichap,bbb,phiwc,betafl,r_0]=fire_ros_new(fuel);
  59 delta_tign=delta_tign_calculation(long,lat,vf,uf,dzdxf,dzdyf,ichap,bbb,phiwc,betafl,r_0);
  60 % Calculates needed variables for rate of fire spread calculation
  61
  62 %%%%%%% First part %%%%%%%
  63 % Set everything inside to time_now and update the tign of the points outside
  64
  65 % Initializing flag matrix A and time of ignition (tign)
  66 % Extending the boundaries, in order to speed up the algorythm
  67 A=[];
  68 C=zeros(n+2,m+2);
  69 % A contains coordinates of the points that were updated during the last
  70 % step
  71 IN_ext=(2)*ones(n+2,m+2);
  72 IN_ext(2:n+1,2:m+1)=IN(:,:,1);
  73
  74 for i=2:n+1
  75     for j=2:m+1
  76         if IN_ext(i,j)==1
  77             if sum(sum(IN_ext(i-1:i+1,j-1:j+1)))<9
  78                 A=[A;[i,j]];
  79             end
  80        end
  81     end
  82 end
  83 tign=ones(n+2,m+2)*1000*time_now;
  84 tign(2:n+1,2:m+1)=IN(:,:,1)*time_now+(1-IN(:,:,1))*1000*time_now;       % and their time of ignition is set to time_now
  85 toc
  86
  87
  88 changed=1;
  89
  90 % The algorithm stops when the matrix converges (tign_old-tign==0) or if
  91 % the amount of iterations
  92 % reaches the max(size()) of the mesh
  93 for istep=1:max(2*size(tign)),
  94     if changed==0,
  95         % The matrix coverged
  96         'first part done'
  97         break
  98     end
  99
 100
 101     tign_last=tign;
 102     time_toc=toc;
 103 str=sprintf('%f -- How long does it take to run step %i',time_toc,istep-1);
 104
 105     % tign_update - updates the time of ignition of the points
 106
 107 [tign,A,C]=tign_update(tign,A,IN_ext,delta_tign,time_now,0);
 108     % tign_update - updates the time of ignition of the points
 109
 110     changed=sum(tign(:)~=tign_last(:));
 111
 112    sprintf('%s \n step %i outside tign changed at %i points \n %f -- norm of the difference',str,istep,changed,norm(tign-tign_last))
 113
 114    end
 115
 116 if changed~=0,
 117
 118    'did not find fixed point'
 119 end
 120
 121 %%%%%%% Second part %%%%%%%
 122
 123 % Set all the points outside to time_now and update the points inside
 124
 125 % Initializing flag matrix A and time of ignition (tign)
 126 % Extending the boundaries, in order to speed up the algorythm
 127 A=[];
 128 C=zeros(n+2,m+2);
 129 for i=2:n+1
 130     for j=2:m+1
 131         if IN_ext(i,j)==0
 132             if sum(sum(IN_ext(i-1:i+1,j-1:j+1)))>0
 133          A=[A;[i,j]];
 134             end
 135        end
 136     end
 137 end
 138
 139 tign_in=zeros(n+2,m+2);
 140 tign_in(2:n+1,2:m+1)=(1-IN(:,:,1)).*tign(2:n+1,2:m+1);
 141 changed=1;
 142
 143 % The algorithm stops when the matrix converges (tign_old-tign==0) or if the amount of iterations
 144 % % reaches the max(size()) of the mesh
 145 for istep=1:max(size(tign)),
 146     if changed==0,
 147                 % The matrix of tign converged
 148                 'printed'
 149                 break
 150     end
 151
 152     tign_last=tign_in;
 153     time_toc=toc;
 154    str= sprintf('%f -- How long does it take to run step %i',time_toc,istep-1);
 155
 156
 157     % tign_update - updates the tign of the points
 158     [tign_in,A,C]=tign_update(tign_in,A,IN_ext,delta_tign,time_now,1);
 159   % hwen it is outside the last parameter is 0, inside 1
 160     changed=sum(tign_in(:)~=tign_last(:));
 161
 162     sprintf('%s \n step %i inside tign changed at %i points \n %f -- norm of the difference',str,istep,changed,norm(tign_in-tign_last))
 163
 164 end
 165 final_tign=zeros(n+2,m+2);
 166 final_tign(2:n+1,2:m+1)=(IN(:,:,1)>0).*tign_in(2:n+1,2:m+1)+(IN(:,:,1)==0).*tign(2:n+1,2:m+1);
 167 result=final_tign(2:n+1,2:m+1);
 168
 169 fid = fopen('output_tign.txt', 'w');
 170     dlmwrite('output_tign.txt', result, 'delimiter', '\t','precision', '%.4f');
 171     fclose(fid);
 172
 173 if changed~=0,
 174     'did not find fixed point inside'
 175    end
 176 end
 177
 178 function [tign,A,C]=tign_update(tign,A,IN,delta_tign,time_now,where)
 179
 180 % Does one iteration of the algorythm and updates the tign of the points of
 181 % the mesh that are next to the previously updated neighbors
 182 %B=A(end,:);
 183 C=zeros(size(tign,1),size(tign,2));
 184 for i=1:size(A,1)
 185     for dx=-1:1
 186         for dy=-1:1
 187                 if where*IN(A(i,1)+dx,A(i,2)+dy)==1
 188                     tign_new=tign(A(i,1),A(i,2))-0.5*(delta_tign(A(i,1)+dx,A(i,2)+dy,dx+2,dy+2)+delta_tign(A(i,1),A(i,2),2-dx,2-dy));
 189                     if (tign(A(i,1)+dx,A(i,2)+dy)<tign_new)&&(tign_new<=time_now)
 190                         % Looking for the max tign, which
 191                         % should be <= than time_now, since the
 192                         % point is inside of the preimeter
 193                    %     if (B(end,1)~=A(i,1)+dx)||(B(end,2)~=A(i,2)+dy)
 194                    %         B=[B;[A(i,1)+dx,A(i,2)+dy]];
 195                             tign(A(i,1)+dx,A(i,2)+dy)=tign_new;
 196                             C(A(i,1)+dx,A(i,2)+dy)=1;
 197                    %    end
 198                     end
 199
 200                 elseif (1-where)*(1-IN(A(i,1)+dx,A(i,2)+dy))==1
 201                     tign_new=tign(A(i,1),A(i,2))+0.5*(delta_tign(A(i,1)+dx,A(i,2)+dy,dx+2,dy+2)+delta_tign(A(i,1),A(i,2),2-dx,2-dy));
 202
 203 if (A(i,1)+dx==2)&&(A(i,2)+dy==2)
 204     display('**********************************************************')
 205 display('i=1,j=1, (2,2) in extended array')
 206 display('we calculate it from the point')
 207 A(i,1)
 208 A(i,2)
 209 display('tign_new that was calculated')
 210 tign_new
 211 display('tign(A(i,1),A(i,2))')
 212 tign(A(i,1),A(i,2))
 213 display('delta_tign(A(i,1)+dx,A(i,2)+dy,dx+2,dy+2) ')
 214 delta_tign(A(i,1)+dx,A(i,2)+dy,dx+2,dy+2)
 215 display('delta_tign(A(i,1),A(i,2),2-dx,2-dy)')
 216 delta_tign(A(i,1),A(i,2),2-dx,2-dy)
 217 display('tign(A(i,1)+dx,A(i,2)+dy')
 218 tign(A(i,1)+dx,A(i,2)+dy)
 219
 220 display('time_now')
 221 time_now
 222 display('check (tign(A(i,1)+dx,A(i,2)+dy)>tign_new)&&(tign_new>=time_now)')
 223 (tign(A(i,1)+dx,A(i,2)+dy)>tign_new)&&(tign_new>=time_now)
 224 end
 225
 226                     if (tign(A(i,1)+dx,A(i,2)+dy)>tign_new)&&(tign_new>=time_now);
 227                         % Looking for the min tign, which
 228                         % should be >= than time_now, since the
 229                         % point is outside of the preimeter
 230                    %     if (B(end,1)~=A(i,1)+dx+1)||(B(end,2)~=A(i,2)+dy)
 231                    %     B=[B;[A(i,1)+dx,A(i,2)+dy]];
 232                         tign(A(i,1)+dx,A(i,2)+dy)=tign_new;
 233                         C(A(i,1)+dx,A(i,2)+dy)=1;
 234                    %     end
 235                     end
 236             end
 237         end
 238     end
 239 end
 240 A=[];
 241 [A(:,1),A(:,2)]=find(C>0);
 242 end
 243
 244 function [ichap,bbb,phiwc,betafl,r_0]=fire_ros_new(fuel,fmc_g)
 245 % ros=fire_ros(fuel,speed,tanphi)
 246 % ros=fire_ros(fuel,speed,tanphi,fmc_g)
 247 % in
 248 %       fuel    fuel description structure
 249 %       speed   wind speed
 250 %       tanphi  slope
 251 %       fmc_g   optional, overrides fuelmc_g from the fuel description
 252 % out
 253 %       ros     rate of spread
 254
 255 % given fuel params
 256
 257 windrf=fuel.windrf;               % WIND REDUCTION FACTOR
 258 fgi=fuel.fgi;                     % INITIAL TOTAL MASS OF SURFACE FUEL (KG/M**2)
 259 fueldepthm=fuel.fueldepthm;       % FUEL DEPTH (M)
 260 savr=fuel.savr;                   % FUEL PARTICLE SURFACE-AREA-TO-VOLUME RATIO, 1/FT
 261 fuelmce=fuel.fuelmce;             % MOISTURE CONTENT OF EXTINCTION
 262 fueldens=fuel.fueldens;           % OVENDRY PARTICLE DENSITY, LB/FT^3
 263 st=fuel.st;                       % FUEL PARTICLE TOTAL MINERAL CONTENT
 264 se=fuel.se;                       % FUEL PARTICLE EFFECTIVE MINERAL CONTENT
 265 weight=fuel.weight;               % WEIGHTING PARAMETER THAT DETERMINES THE SLOPE OF THE MASS LOSS CURVE
 266 fci_d=fuel.fci_d;                 % INITIAL DRY MASS OF CANOPY FUEL
 267 fct=fuel.fct;                     % BURN OUT TIME FOR CANOPY FUEL, AFTER DRY (S)
 268 ichap=fuel.ichap;                 % 1 if chaparral, 0 if not
 269 fci=fuel.fci;                     % INITIAL TOTAL MASS OF CANOPY FUEL
 270 fcbr=fuel.fcbr;                   % FUEL CANOPY BURN RATE (KG/M**2/S)
 271 hfgl=fuel.hfgl;                   % SURFACE FIRE HEAT FLUX THRESHOLD TO IGNITE CANOPY (W/m^2)
 272 cmbcnst=fuel.cmbcnst;             % JOULES PER KG OF DRY FUEL
 273 fuelheat=fuel.fuelheat;           % FUEL PARTICLE LOW HEAT CONTENT, BTU/LB
 274 fuelmc_g=fuel.fuelmc_g;           % FUEL PARTICLE (SURFACE) MOISTURE CONTENT, jm: 1 by weight?
 275 fuelmc_c=fuel.fuelmc_c;           % FUEL PARTICLE (CANOPY) MOISTURE CONTENT, 1
 276
 277 if exist('fmc_g','var') % override moisture content by given
 278     fuelmc_g = fmc_g;
 279 end
 280
 281 % computations from CAWFE code: wf2_janice/fire_startup.m4
 282
 283 bmst     = fuelmc_g/(1+fuelmc_g);          % jm: 1
 284 fuelheat = cmbcnst * 4.30e-04;             % convert J/kg to BTU/lb
 285 fci      = (1.+fuelmc_c)*fci_d;
 286 fuelloadm= (1.-bmst) * fgi;                % fuelload without moisture
 287                                            % jm: 1.-bmst = 1/(1+fuelmc_g) so fgi includes moisture?
 288 fuelload = fuelloadm * (.3048)^2 * 2.205;  % to lb/ft^2
 289 fueldepth= fueldepthm/0.3048;              % to ft
 290 betafl   = fuelload/(fueldepth * fueldens);% packing ratio  jm: lb/ft^2/(ft * lb*ft^3) = 1
 291 betaop   = 3.348 * savr^(-0.8189);         % optimum packing ratio jm: units??
 292 qig      = 250. + 1116.*fuelmc_g;          % heat of preignition, btu/lb
 293 epsilon  = exp(-138./savr );               % effective heating number
 294 rhob     = fuelload/fueldepth;             % ovendry bulk density, lb/ft^3
 295 c        = 7.47 * exp(-0.133 * savr^0.55); % const in wind coef
 296 bbb      = 0.02526 * savr^0.54;            % const in wind coef
 297 c        = c * windrf^bbb;                 % jm: wind reduction from 20ft per Baughman&Albini(1980)
 298 e        = 0.715 * exp( -3.59e-4 * savr);  % const in wind coef
 299 phiwc    = c * (betafl/betaop)^(-e);
 300 rtemp2   = savr^1.5;
 301 gammax   = rtemp2/(495. + 0.0594*rtemp2);  % maximum rxn vel, 1/min
 302 a        = 1./(4.774 * savr^0.1 - 7.27);   % coef for optimum rxn vel
 303 ratio    = betafl/betaop;
 304 gamma    = gammax*(ratio^a)*exp(a*(1.-ratio)); % optimum rxn vel, 1/min
 305 wn       = fuelload/(1 + st);              % net fuel loading, lb/ft^2
 306 rtemp1   = fuelmc_g/fuelmce;
 307 etam     = 1.-2.59*rtemp1 +5.11*rtemp1^2 -3.52*rtemp1^3;  % moist damp coef
 308 etas     = 0.174* se^(-0.19);              % mineral damping coef
 309 ir       = gamma * wn * fuelheat * etam * etas; % rxn intensity,btu/ft^2 min
 310 irm      = ir * 1055./( 0.3048^2 * 60.) * 1.e-6;% for mw/m^2 (set but not used)
 311 xifr     = exp( (0.792 + 0.681*savr^0.5)...
 312             * (betafl+0.1)) /(192. + 0.2595*savr); % propagating flux ratio
 313 %        ... r_0 is the spread rate for a fire on flat ground with no wind.
 314 r_0      = ir*xifr/(rhob * epsilon *qig);  % default spread rate in ft/min
 315
 316 % computations from CAWFE code: wf2_janice/fire_ros.m4
 317
 318
 319 end
 320
 321 function delta_tign=delta_tign_calculation(long,lat,vf,uf,dzdxf,dzdyf,ichap,bbb,phiwc,betafl,r_0)
 322     %Extend the boundaries to speed up the algorithm, the values of the
 323     %extended boundaries would be set to zeros and are never used in the
 324     %code
 325         result=1;
 326     'hello'
 327     delta_tign=zeros(size(long,1)+2,size(long,2)+2,3,3);
 328     rate_of_spread=zeros(size(long,1)+2,size(long,2)+2,3,3);
 329
 330     long2=zeros(size(long,1)+2,size(long,2)+2);
 331     long2(2:size(long,1)+1,2:size(long,2)+1)=long;
 332     long=long2;
 333
 334     lat2=zeros(size(lat,1)+2,size(lat,2)+2);
 335     lat2(2:size(lat,1)+1,2:size(lat,2)+1)=lat;
 336     lat=lat2;
 337
 338     vf2=zeros(size(vf,1)+2,size(vf,2)+2);
 339     vf2(2:size(vf,1)+1,2:size(vf,2)+1)=vf;
 340     vf=vf2;
 341
 342     uf2=zeros(size(uf,1)+2,size(uf,2)+2);
 343     uf2(2:size(uf,1)+1,2:size(uf,2)+1)=uf;
 344     uf=uf2;
 345
 346     dzdxf2=zeros(size(dzdxf,1)+2,size(dzdxf,2)+2);
 347     dzdxf2(2:size(dzdxf,1)+1,2:size(dzdxf,2)+1)=dzdxf;
 348     dzdxf=dzdxf2;
 349
 350     dzdyf2=zeros(size(dzdyf,1)+2,size(dzdyf,2)+2);
 351     dzdyf2(2:size(dzdyf,1)+1,2:size(dzdyf,2)+1)=dzdyf;
 352     dzdyf=dzdyf2;
 353
 354 for i=2:size(long,1)-1
 355     for j=2:size(long,2)-1
 356         for a=-1:1
 357             for b=-1:1
 358                 wind=0.5*((long(i,j,1)-long(i+a,j+b,1))*vf(i,j,1)+  ...
 359                       (lat(i,j,1)-lat(i+a,j+b,1))*uf(i,j,1));
 360                 angle=0.5*((long(i,j,1)-long(i+a,j+b,1))*dzdxf(i,j,1)+  ...
 361                        (lat(i,j,1)-lat(i+a,j+b,1))*dzdyf(i,j,1));
 362                 if ~ichap,
 363                     %       ... if wind is 0 or into fireline, phiw = 0, &this reduces to backing ros.
 364                     spdms = max(wind,0.);
 365                     umidm = min(spdms,30.);                    % max input wind spd is 30 m/s   !param!
 366                     umid = umidm * 196.850;                    % m/s to ft/min
 367                     %  eqn.: phiw = c * umid**bbb * (betafl/betaop)**(-e) ! wind coef
 368                     phiw = umid^bbb * phiwc;                   % wind coef
 369                     phis = 5.275 * betafl^(-0.3) *max(0,angle)^2;   % slope factor
 370                     ros = r_0*(1. + phiw + phis)  * .00508; % spread rate, m/s
 371                 else  % chapparal
 372                     %        .... spread rate has no dependency on fuel character, only windspeed.
 373                     spdms = max(wind,0.);
 374                     ros = max(.03333,1.2974 * spdms^1.41);       % spread rate, m/s
 375                 end
 376                 rate_of_spread(i,j,a+2,b+2)=min(ros,6);
 377                 % DEscribe the coefficient below
 378                 delta_tign(i,j,a+2,b+2)=sqrt((long(i+a,j+b,1)-long(i,j,1))^2+    ...
 379                           (lat(i+a,j+b,1)-lat(i,j,1))^2)/rate_of_spread(i,j,a+2,b+2);
 380             end
 381         end
 382
 383     end
 384
 385 end
 386 end
 387
 388 function [tign,A]=point_update(i,j,tign,delta_tign,inside)
 389
 390     for dx=-1:+1
 391         for dy=-1:+1
 392             % loop over all neighbors
 393                 if (A(i+dx,j+dy)==1) % the neighbor was already updated
 394                 % All the vectors are split in half-intervals
 395                 % to get better calculation
 396             %%% Make a picture of what is happening %%%
 397
 398             tign_new=tign(i+dx,j+dy)-delta_tign(i,j,dx+2,dy+2)-delta_tign(i+dx,j+dy,2-dx,2-dy);
 399
 400                 if (tign(i,j)<tign_new)&&(tign_new<=time_now);
 401                     % Looking for the max tign, which
 402                     % should be <= than time_now, since the
 403                     % point is inside of the preimeter
 404                     tign(i,j)=tign_new;
 405                     A(i,j)=1;
 406                 end
 407             end
 408         end
 409     end
 410 end
 411
 412 % function result=print_matrix(tign,fid)
 413 %
 414 %
 415 % fprintf(fid,'%s\n','j=1740:1744, i=2750:2754');
 416 %
 417 % for ii = 2750:2754
 418 %     fprintf(fid,'%g\t',tign(ii,1740:1744));
 419 %     fprintf(fid,'\n');
 420 % end
 421 % fprintf(fid,'%s\n','i=100');
 422 %
 423 % for ii = 98:102
 424 %     fprintf(fid,'%g\t',tign(ii,98:102));
 425 %     fprintf(fid,'\n');
 426 % end
 427 % fprintf(fid,'%s\n','i=1000');
 428 %
 429 % for ii = 998:1002
 430 %     fprintf(fid,'%g\t',tign(ii,998:1002));
 431 %     fprintf(fid,'\n');
 432 % end
 433 % result=0;
 434 % end
 435
 436 % for i=2:size(tign,1)-1
 437 %     for j=2:size(tign,2)-1
 438 %         % sum_A is needed to know what is the amount of points that surrounds (i,j)
 439 %         sum_A=sum(sum(A(i-1:i+1,j-1:j+1)));
 440 %
 441 %         if (sum_A~=0)
 442 %
 443 %             % sum_A>0 then at least on eneighbor was previously updated and its
 444 %             % tign can be used to update the tign of the point (i,j)
 445 %
 446 %             % I subtract 1 in all IN, long, lat, uf, vf, dzdxf, dzdyf
 447 %             % matrices, because their boundaries were not updated, unlike
 448 %             %%% tign (do the loop 1:n, 1:m and change the indexes in tign) %%%
 449 %
 450 %
 451 %             for dx=-1:1
 452 %                 for dy=-1:1
 453 %                     % loop over all neighbors
 454 %                     if (A(i+dx,j+dy)==1) % the neighbor was already updated
 455 %                         % All the vectors are split in half-intervals
 456 %                         % to get better calculation
 457 %                         %%% Make a picture of what is happening %%%
 458 %                         % I do multiplication by 0.5 here 0.5-(IN(i-1,j-1)>0
 459 %
 460 %                         tign_new=tign(i+dx,j+dy)+(0.5-(IN(i-1,j-1)>0))*(delta_tign(i,j,dx+2,dy+2)+delta_tign(i+dx,j+dy,2-dx,2-dy));
 461 %
 462 %                         if (IN(i-1,j-1)>0)
 463 %
 464 %                             if (tign(i,j)<tign_new)&&(tign_new<=time_now)
 465 %                                 % Looking for the max tign, which
 466 %                                 % should be <= than time_now, since the
 467 %                                 % point is inside of the preimeter
 468 %                                 tign(i,j)=tign_new;
 469 %                                 A(i,j)=1;
 470 %                             end
 471 %                         elseif (tign(i,j)>tign_new)&&(tign_new>=time_now);
 472 %                                 % Looking for the min tign, which
 473 %                                 % should be >= than time_now, since the
 474 %                                 % point is outside of the preimeter
 475 %                                 tign(i,j)=tign_new;
 476 %                                 A(i,j)=1;
 477 %                             end
 478 %                         end
 479 %                     end
 480 %                 end
 481 %             end
 482 %         end
 483 %     end
 484
 485
 486
 487 % Old Code -- tign_update
 488
 489
 490 %                             % wind1 = vect.*(uf,vf)
 491 %                             wind1=0.5*((long(a-1,b-1,1)-long(i-1,j-1, 1))*vf(i-1,j-1,1)+  ...
 492 %                                      (lat(a-1,b-1,1)-lat(i-1,j-1,1))*uf(i-1,j-1,1));
 493 %                             angle1=0.5*((long(a-1,b-1,1)-long(i-1,j-1,1))*dzdxf(i-1,j-1,1)+  ...
 494 %                                      (lat(a-1,b-1,1)-lat(i-1,j-1,1))*dzdyf(i-1,j-1,1));
 495 %
 496 %                             % This is needed to calculate the ros(i,j)
 497 %                              if ~ichap,
 498 %                                 %       ... if wind is 0 or into fireline, phiw = 0, &this reduces to backing ros.
 499 %                                 spdms = max(wind1,0.);
 500 %                                 umidm = min(spdms,30.);                    % max input wind spd is 30 m/s   !param!
 501 %                                 umid = umidm * 196.850;                    % m/s to ft/min
 502 %                                 %  eqn.: phiw = c * umid**bbb * (betafl/betaop)**(-e) ! wind coef
 503 %                                 phiw = umid^bbb * phiwc;                   % wind coef
 504 %                                 phis = 5.275 * betafl^(-0.3) *max(0,angle1)^2;   % slope factor
 505 %                                 ros = r_0*(1. + phiw + phis)  * .00508; % spread rate, m/s
 506 %                             else  % chapparal
 507 %                                 %        .... spread rate has no dependency on fuel character, only windspeed.
 508 %                                 spdms = max(wind1,0.);
 509 %                                 ros = max(.03333,1.2974 * spdms^1.41);       % spread rate, m/s
 510 %                             end
 511 %                             ros=min(ros,6);
 512 %
 513 %                             % tign_new=tign(a,b)-delta(t);
 514 %                               tign_new1=tign(a,b)-0.5*sqrt((long(a-1,b-1,1)-long(i-1,j-1,1))^2+   ...
 515 %                                      (lat(a-1,b-1,1)-lat(i-1,j-1,1))^2)/ros;                   ...
 516 %
 517 %                             % Same calculation for the second half of the
 518 %                             % interval
 519 %
 520 %                             wind2=0.5*((long(a-1,b-1,1)-long(i-1,j-1, 1))*vf(a-1,b-1,1)+  ...
 521 %                                      (lat(a-1,b-1,1)-lat(i-1,j-1,1))*uf(a-1,b-1,1));
 522 %                             angle2=0.5*((long(a-1,b-1,1)-long(i-1,j-1,1))*dzdxf(a-1,b-1,1)+  ...
 523 %                                      (lat(a-1,b-1,1)-lat(i-1,j-1,1))*dzdyf(a-1,b-1,1));
 524 %                             if ~ichap,
 525 %                                 %       ... if wind is 0 or into fireline, phiw = 0, &this reduces to backing ros.
 526 %                                 spdms = max(wind2,0.);
 527 %                                 umidm = min(spdms,30.);                    % max input wind spd is 30 m/s   !param!
 528 %                                 umid = umidm * 196.850;                    % m/s to ft/min
 529 %                                 %  eqn.: phiw = c * umid**bbb * (betafl/betaop)**(-e) ! wind coef
 530 %                                 phiw = umid^bbb * phiwc;                   % wind coef
 531 %                                 phis = 5.275 * betafl^(-0.3) *max(0,angle2)^2;   % slope factor
 532 %                                 ros = r_0*(1. + phiw + phis)  * .00508; % spread rate, m/s
 533 %                             else  % chapparal
 534 %                                 %        .... spread rate has no dependency on fuel character, only windspeed.
 535 %                                 spdms = max(wind2,0.);
 536 %                                 ros = max(.03333,1.2974 * spdms^1.41);       % spread rate, m/s
 537 %                             end
 538 %                             ros=min(ros,6);
 539 %                               tign_new2=-0.5*sqrt((long(a-1,b-1,1)-long(i-1,j-1,1))^2+   ...
 540 %                                      (lat(a-1,b-1,1)-lat(i-1,j-1,1))^2)/ros;
 541 %                             tign_new=tign_new1+tign_new2;
 542
 543
 544
 545 % OUTSIDE
 546
 547
 548
 549
 550 %                             wind1=0.5*((long(i-1,j-1,1)-long(a-1,b-1,1))*vf(i-1,j-1,1)+  ...
 551 %                                      (lat(i-1,j-1,1)-lat(a-1,b-1,1))*uf(i-1,j-1,1));
 552 %                             angle1=0.5*((long(i-1,j-1,1)-long(a-1,b-1,1))*dzdxf(i-1,j-1,1)+  ...
 553 %                                      (lat(i-1,j-1,1)-lat(a-1,b-1,1))*dzdyf(i-1,j-1,1));
 554 %                             if ~ichap,
 555 %                                 %       ... if wind is 0 or into fireline, phiw = 0, &this reduces to backing ros.
 556 %                                 spdms = max(wind1,0.);
 557 %                                 umidm = min(spdms,30.);                    % max input wind spd is 30 m/s   !param!
 558 %                                 umid = umidm * 196.850;                    % m/s to ft/min
 559 %                                 %  eqn.: phiw = c * umid**bbb * (betafl/betaop)**(-e) ! wind coef
 560 %                                 phiw = umid^bbb * phiwc;                   % wind coef
 561 %                                 phis = 5.275 * betafl^(-0.3) *max(0,angle1)^2;   % slope factor
 562 %                                 ros = r_0*(1. + phiw + phis)  * .00508; % spread rate, m/s
 563 %                             else  % chapparal
 564 %                                 %        .... spread rate has no dependency on fuel character, only windspeed.
 565 %                                 spdms = max(wind1,0.);
 566 %                                 ros = max(.03333,1.2974 * spdms^1.41);       % spread rate, m/s
 567 %                             end
 568 %                             ros=min(ros,6);
 569 %                             % tign_new=tign(a,b)-delta(t);
 570 %                             tign_new1=tign(a,b)+0.5*sqrt((long(a-1,b-1,1)-long(i-1,j-1,1))^2+    ...
 571 %                                      (lat(a-1,b-1,1)-lat(i-1,j-1,1))^2)/r
 572 %                                      os;s
 573 %
 574 %                             wind2=0.5*((long(i-1,j-1,1)-long(a-1,b-1,1))*vf(a-1,b-1,1)+  ...
 575 %                                      (lat(i-1,j-1,1)-lat(a-1,b-1,1))*uf(a-1,b-1,1));
 576 %                             angle2=0.5*((long(i-1,j-1,1)-long(a-1,b-1,1))*dzdxf(a-1,b-1,1)+  ...
 577 %                                      (lat(i-1,j-1,1)-lat(a-1,b-1,1))*dzdyf(a-1,b-1,1));
 578 %                             if ~ichap,
 579 %                                 %       ... if wind is 0 or into fireline, phiw = 0, &this reduces to backing ros.
 580 %                                 spdms = max(wind2,0.);
 581 %                                 umidm = min(spdms,30.);                    % max input wind spd is 30 m/s   !param!
 582 %                                 umid = umidm * 196.850;                    % m/s to ft/min
 583 %                                 %  eqn.: phiw = c * umid**bbb * (betafl/betaop)**(-e) ! wind coef
 584 %                                 phiw = umid^bbb * phiwc;                   % wind coef
 585 %                                 phis = 5.275 * betafl^(-0.3) *max(0,angle2)^2;   % slope factor
 586 %                                 ros = r_0*(1. + phiw + phis)  * .00508; % spread rate, m/s
 587 %                             else  % chapparal
 588 %                                 %        .... spread rate has no dependency on fuel character, only windspeed.
 589 %                                 spdms = max(wind2,0.);
 590 %                                 ros = max(.03333,1.2974 * spdms^1.41);       % spread rate, m/s
 591 %                             end
 592 %                             ros=min(ros,6);
 593 %                             tign_new2=0.5*sqrt((long(a-1,b-1,1)-long(i-1,j-1,1))^2+    ...
 594 %                                      (lat(a-1,b-1,1)-lat(i-1,j-1,1))^2)/ros;
 595 %                             tign_new=tign_new1+tign_new2;                    ...
 596 %                             % Here the direction of the vector is
 597 %                             % opposite, since fire is going from the
 598 %                             % inside point towards the point that was
 599 %                             % already computed
 600
 601
 602
 603
 604
 605
 606
 607
 608
 609
 610
 611
 612
 613
 614
 615
 616