%这是主程序，调用其他子程序
clear all;
% 下面的kmax表示最大度，p表示分布。都是把BA模型中计算的结果先保存下来，然后复制到这里。（每次得到结果不一样。因为BA模型中有随机因素）
kmax =82;
p =[0 0 0.3940 0.1950 0.1180 0.0750 0.0540 0.0260 0.0250 0.0180 0.0190 0.0110 0.0080 0.0030 0.0080 0.0060 0.0030 0.0030 0.0030 0 0.0040 0.0010 0.0020 0.0030 0 0.0020 0 0.0020 0.0010 0.0010 0 0 0.0020 0 0 0 0.0010 0.0010 0.0010 0 0 0.0020 0.0020 0 0 0 0 0 0 0 0.0020 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0020 0 0 0 0 0 0 0 0 0 0.0010 0 0 0 0 0 0 0.0010];
%度分布可以通过BA模型生成，我的下载资源里面有BA模型。
lambda=0;% 传染率
%循环是不断的增加lambda,看lambda何值时会出现疾病爆发
for i=1:100
[t,I] = ode45('targetimmunity_c',[0,20],rand(kmax,1)*0.01,[],lambda);
%调用其他子程序，[0,20]是时间初值，rand(kmax,1)*0.01是赋被感染者随机初值，[]是格式要求.
size=length(I(:,1));
B=0;%B是为了计算每个lambda对应的最终被感染人数的比例。
for k=1:kmax
B=B+I(size,k)*p(k);
end
%T 和 x的作用是为了把结果以向量的形式输出。
T(i)=B;
x(i)=lambda;
lambda=lambda+0.01;
end
plot(x,T,'.-r');
hold on;
