close all clear clc filename = ['wea.csv']; t1 = dataset('file',filename,'delimiter',',','format','%s %d %d %d %f %f'); t1.Date = datenum(t1.Date_Time,'yyyy-mm-dd'); t2 = t1(ismember(t1.Month,[5:10]),:); t2.ddm = ones(length(t2),1); for i = 32:184 t2.ddm(i) = sum(t2.MeanTemp__C_(i-11:i-1))/11; end t3 = t2(ismember(t2.Month,[6,7,8,9,10]),:); tt = 1:length(t3); Tspan = [1,160]; y0 = [105,108,108,0,1,6,0,0,0]; %% two cases of SWMP [time0, Y0] = ode45(@(t,y) WeatherCases(t,y,tt,t3.MeanTemp__C_,t3.ddm,t3.TotalPrecip_mm_,t3.TotalPrecip_mm_,9,1),Tspan,y0); % Solve ODE [time, Y] = ode45(@(t,y) WeatherCases(t,y,tt,t3.MeanTemp__C_,t3.ddm,t3.TotalPrecip_mm_,t3.TotalPrecip_mm_,9,2),Tspan,y0); % Solve ODE figure plot(time0,Y0(:,2)+Y0(:,3)+Y0(:,4)+Y0(:,5),'color','b','LineWidth',1.5) hold on plot(time,Y(:,2)+Y(:,3)+Y(:,4)+Y(:,5),'color','g','LineWidth',1.5) xlabel('t') ylabel('Total adult mosquito population') legend('Case 1','Case 2') hold off figure [AX0,H10,H20] = plotyy(time0,Y0(:,5),time0,Y0(:,7),'plot'); set(H10,'LineStyle','--','LineWidth',1.5,'color','b') set(H20,'LineStyle','--','LineWidth',1.5,'color','g') hold(AX0(1)) H1 = plot(AX0(1),time,Y(:,5),'LineStyle','-','LineWidth',1.5,'color','b'); hold(AX0(2)) H2 = plot(AX0(2),time,Y(:,7),'LineStyle','-','LineWidth',1.5,'color','g'); set(get(AX0(1),'Ylabel'),'String','Infectious female mosquito population','color','b') set(get(AX0(2),'Ylabel'),'String','Infectiou bird population','color','g') xlabel('t') set(AX0(1),'ycolor','b') set(AX0(2),'ycolor','g') hold off legend([H10;H20;H1;H2],'Case 1','Case 1','Case 2','Case 2') %% impacts of rho [time4, Y4] = ode45(@(t,y) Weather(t,y,tt,t3.MeanTemp__C_,t3.ddm,t3.TotalPrecip_mm_,t3.TotalPrecip_mm_,9,0.0001),Tspan,y0); % Solve ODE [time5, Y5] = ode45(@(t,y) Weather(t,y,tt,t3.MeanTemp__C_,t3.ddm,t3.TotalPrecip_mm_,t3.TotalPrecip_mm_,29,0.0001),Tspan,y0); % Solve ODE figure plot(time4,Y4(:,2)+Y4(:,3)+Y4(:,4)+Y4(:,5),'color','b','LineWidth',1.5) hold on plot(time5,Y5(:,2)+Y5(:,3)+Y5(:,4)+Y5(:,5),'color','g','LineWidth',1.5) xlabel('t') ylabel('Total adult mosquito population') legend('\rho=1','\rho=29') hold off figure [AX4,H14,H24] = plotyy(time4,Y4(:,5),time4,Y4(:,7),'plot'); set(H14,'LineStyle','--','LineWidth',1.5,'color','b') set(H24,'LineStyle','--','LineWidth',1.5,'color','g') hold(AX4(1)) H15 = plot(AX4(1),time5,Y5(:,5),'LineStyle','-','LineWidth',1.5,'color','b'); hold(AX4(2)) H25 = plot(AX4(2),time5,Y5(:,7),'LineStyle','-','LineWidth',1.5,'color','g'); set(get(AX4(1),'Ylabel'),'String','Infectious female mosquito population','color','b') set(get(AX4(2),'Ylabel'),'String','Infectiou bird population','color','g') xlabel('t') set(AX4(1),'ycolor','b') set(AX4(2),'ycolor','g') hold off legend([H14;H24;H15;H25],'\rho=1','\rho=1','\rho=29','\rho=29') %% the impact of kappa [time1, Y1] = ode45(@(t,y) WeatherKappa(t,y,tt,t3.MeanTemp__C_,t3.ddm,t3.TotalPrecip_mm_,t3.TotalPrecip_mm_,9,0),Tspan,y0); % Solve ODE [time2, Y2] = ode45(@(t,y) WeatherKappa(t,y,tt,t3.MeanTemp__C_,t3.ddm,t3.TotalPrecip_mm_,t3.TotalPrecip_mm_,9,0.001),Tspan,y0); % Solve ODE figure plot(time1,Y1(:,2)+Y1(:,3)+Y1(:,4)+Y1(:,5),'color','b','LineWidth',1.5) hold on plot(time2,Y2(:,2)+Y2(:,3)+Y2(:,4)+Y2(:,5),'color','g','LineWidth',1.5) xlabel('t') ylabel('Total adult mosquito population') h = legend('$\overline{\kappa}=0$','$\overline \kappa=0.001$' ); set(h,'Interpreter','latex') figure [AX1,H11,H12] = plotyy(time1,Y1(:,5),time1,Y1(:,7),'plot'); set(H11,'LineStyle','--','LineWidth',1.5,'color','b') set(H12,'LineStyle','--','LineWidth',1.5,'color','g') hold(AX1(1)) H21 = plot(AX1(1),time2,Y2(:,5),'LineStyle','-','LineWidth',1.5,'color','b'); hold(AX1(2)) H22 = plot(AX1(2),time2,Y2(:,7),'LineStyle','-','LineWidth',1.5,'color','g'); set(get(AX1(1),'Ylabel'),'String','Infectious female mosquito population','color','b') set(get(AX1(2),'Ylabel'),'String','Infectiou bird population','color','g') xlabel('t') set(AX1(2),'ycolor','g') set(AX1(1),'ycolor','b') hold off hh = legend([H11;H12;H21;H22],'$\overline{\kappa}=0$','$\overline{\kappa}=0$','$\overline \kappa=0.001$','$\overline{\kappa}=0.001$'); set(hh,'Interpreter','latex') %% The impact of rain % Different P in July kp = 0.0007; t3.supp1 = zeros(length(t3),1); t3.supp1 (t3.Month == 7) = 30; t3.supp2 = zeros(length(t3),1); t3.supp2 (t3.Month == 7) = 60; P1 = t3.TotalPrecip_mm_; P2 = t3.TotalPrecip_mm_ + t3.supp1; P3 = t3.TotalPrecip_mm_ + t3.supp2; P = [P1,P2,P3]; C = {'r','b','g','m','y',[.5 .6 .7],[.8 .2 .6]}; % Cell array of colros. figure for i = 1:3 PP = P(:,i); [time, Y] = ode45(@(t,y) Weather(t,y,tt,t3.MeanTemp__C_,t3.ddm,PP,t3.TotalPrecip_mm_,9,kp),Tspan,y0); % Solve ODE plot(time,Y(:,2)+Y(:,3)+Y(:,4)+Y(:,5),'color',C{i},'LineWidth',2) hold on xlabel('t') ylabel('Total adult mosquito population') legend( 'Normal precipitation','Heavy precipitation (30mm more)','Heavier precipitation (60mm more)') end figure for i = 1:3 PP = P(:,i); [time, Y] = ode45(@(t,y) Weather(t,y,tt,t3.MeanTemp__C_,t3.ddm,PP,t3.TotalPrecip_mm_,9,kp),Tspan,y0); % Solve ODE plot(time,Y(:,5),'color',C{i},'LineWidth',1.5) hold on xlabel('t') ylabel('Infectious female mosquito population') legend( 'Normal precipitation','Heavy precipitation (30mm more)','Heavier precipitation (60mm more)') end figure for i = 1:3 PP = P(:,i); [time, Y] = ode45(@(t,y) Weather(t,y,tt,t3.MeanTemp__C_,t3.ddm,PP,t3.TotalPrecip_mm_,9,kp),Tspan,y0); % Solve ODE plot(time,Y(:,7),'color',C{i},'LineWidth',1.5) hold on xlabel('t') ylabel('Infectiou bird population') legend( 'Normal precipitation','Heavy precipitation (30mm more)','Heavier precipitation (60mm more)','all more rain') end % Different P in Sept. kp = 0.0007; t3.supp1 = zeros(length(t3),1); t3.supp1 (t3.Month == 9) = 30; t3.supp2 = zeros(length(t3),1); t3.supp2 (t3.Month == 9) = 60; t3.supp3 = 30*ones(length(t3),1); P1 = t3.TotalPrecip_mm_; P2 = t3.TotalPrecip_mm_ + t3.supp1; P3 = t3.TotalPrecip_mm_ + t3.supp2; P = [P1,P2,P3]; C = {'r','b','g','m','y',[.5 .6 .7],[.8 .2 .6]}; % Cell array of colros. figure for i = 1:3 PP = P(:,i); [time, Y] = ode45(@(t,y) Weather(t,y,tt,t3.MeanTemp__C_,t3.ddm,PP,t3.TotalPrecip_mm_,9,kp),Tspan,y0); % Solve ODE plot(time,Y(:,2)+Y(:,3)+Y(:,4)+Y(:,5),'color',C{i},'LineWidth',2) hold on xlabel('t') ylabel('Total adult mosquito population') legend( 'Normal precipitation','Heavy precipitation (30mm more)','Heavier precipitation (60mm more)') end figure for i = 1:3 PP = P(:,i); [time, Y] = ode45(@(t,y) Weather(t,y,tt,t3.MeanTemp__C_,t3.ddm,PP,t3.TotalPrecip_mm_,9,kp),Tspan,y0); % Solve ODE plot(time,Y(:,5),'color',C{i},'LineWidth',1.5) hold on xlabel('t') ylabel('Infectious female mosquito population') legend( 'Normal precipitation','Heavy precipitation (30mm more)','Heavier precipitation (60mm more)') end figure for i = 1:3 PP = P(:,i); [time, Y] = ode45(@(t,y) Weather(t,y,tt,t3.MeanTemp__C_,t3.ddm,PP,t3.TotalPrecip_mm_,9,kp),Tspan,y0); % Solve ODE plot(time,Y(:,7),'color',C{i},'LineWidth',1.5) hold on xlabel('t') ylabel('Infectiou bird population') legend( 'Normal precipitation','Heavy precipitation (30mm more)','Heavier precipitation (60mm more)','all more rain') end