元胞自动机之森林火灾升级版MATLAB代码

元胞自动机中经典森林火灾模型（(1)正在燃烧的树变成空格位； (2)如果绿树格位的最近邻居中有一个树在燃烧，则它变成正在燃烧的树； (3)在空格位，树以概率p生长； (4)在最近的邻居中没有正在燃烧的树的情况下，树在每一步，以概率f(闪电)变为正在燃烧的树。）的拓展：初始状态只有一个着火点，不存在闪电袭击导致着火，从这一个着火点开始蔓延，蔓延的概率达到多大的时候，整个森林可能会被烧光？ MATLAB代码 m文件，创新点是考虑了风向的影响，受到表扬

%矩阵中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;