%矩陣中0代表空地，1代表綠樹，2代表燃燒
%網格圖中白色代表空地，綠色代表綠樹，紅色代表燃燒
%rule?1正在燃燒的樹變成空格位；
%rule?2如果綠樹格位的最近鄰居中有一個樹在燃燒，則它以蔓延概率變成正在燃燒的樹；
%rule?3在空格位，樹以概率p生長；
%%?準備和初始化
close?all;
clc;
clear;
figure;
p=0.001;%s生長概率
f=0.01;%著火概率
M=50;%森林大小
T=100;%運行步數
%蔓延概率
p0=0.3;%風力影響
a=pi/4;%風向夾角
q=1;%計算概率時的臨時變量初值

for?p1=0.75:0.05:01%改變蔓延概率p1大小
????p1
????%%?概率矩陣
????%計算八鄰域內著火的樹木在風力影響下使中心綠樹著火的概率，風力影響通過沿八方向向量分解表達
????P=[];
????g=[cos（a）;sin（a）];
????for?i=1:3;
????????for?j=1:3;
????????????P（ij）=[i-2j-2]*g*p0/sqrt（（i-2）^2+（j-2）^2）+p1;
????????end
????end
????P（22）=0;
????%%?每次實驗賦初值
????flag=0;
????for?l=1:10%對每個蔓延概率，進行十次實驗，記錄“燒光”的次數
????????set（gcf‘DoubleBuffer‘‘on‘）;
????????I=randi（2M）-1;%Inside森林矩陣
????????I（randi（M^2））=2;
????????D=ones（M+2M+23）;%Draw畫RGB森林矩陣
????????D（2:M+12:M+11）=（I~=1）;%R
????????D（2:M+12:M+12）=（I~=2）;%G
????????D（2:M+12:M+13）=（I==0）;%B
????????imshow（D‘InitialMagnification‘‘fit‘）;%使圖像以合適大小顯示
????????t=0;
????????title（[‘p=‘num2str（p1）‘?n=‘num2str（l）‘?T=‘num2str（t）]）;%%時間記錄
????????pause（0.2）

????????for?t=1:T
????????????NI=I;%New?Inside新森林矩陣
????????????NI（I==2）=0;%for?rule?1
????????????NI=NI+（I==0）.*（rand（M）????????????%%?for?rule?2
????????????O=zeros（M+2）;%Outside加邊森林矩陣
????????????O（2:M+12:M+1）=I;
????????????NO=O;%New?Outside新加邊森林矩陣
????????????NO（2:M+12:M+1）=NI;
??