% THE ARMAGEDDON EQUATIONS, Mark VI
%
% simulates the stability/instability of the application of arms control methods
% on a state similar to Pakistan.
%
% PHY 256, Winter 2009, Juliette Zerick (jnzerick@math.ucdavis.edu)

clear all;
colordef black;
format long;

% ENABLE/DISABLE OPTIONAL THINGS
	
ENABLE_BRAKES = 1; % quit iteration if HUGE_VAL is reached
	HUGE_VAL = 50;

% SCALING TERMS	
	
ECON_SANC = 0; % 0: no economic sanctions
MIL_SANC = 0; % 0: no military sanctions
	
% ITERATIONS

numiter=150; % maximum number of iterations

% PLOT SUBSET

PLOT_SUBSET = 0; % 0: switch off; 1: switch on
	pss_start = 1;
	pss_end = 14;

% INITIAL VALUES

	ivscaling = .1;

	G_E(1) = 5		* ivscaling;
	G_S(1) = 7		* ivscaling;
	B_E(1) = 5		* ivscaling;
	B_S(1) = 7		* ivscaling;
	P(1) = 0		* ivscaling;
	M(1) = 6		* ivscaling;
	T(1) = 1		* ivscaling; 

% PARAMETERS

	pscaling = 0.1;

	% G_E
	u_gegs = 3	*pscaling; % support of state
	u_get = 2	*pscaling; % effects of stability
	u_gem = 2	*pscaling; % cost of the military
	
	u_gece = 4	*pscaling; % competition with B_E
	u_gek = 5	*pscaling*ECON_SANC; % economic sanctions

	% G_S
	u_gsge = 3	*pscaling; % strong economy -> strong government
	u_gsm = 4	*pscaling; % security from military/police
	
	u_gsp = 2	*pscaling; % support by the people
	u_gst = 2	*pscaling; % effects of stability

	u_gscs = 3	*pscaling; % competition with B_S
	
	% B_E
	u_bebs = 3	*pscaling; % support of B_S

	u_bet = 2	*pscaling; % effects of stability
	
	u_becs = 2	*pscaling; % crackdown by G_S
	u_bece = 4	*pscaling; % competition with G_E
	u_bek = 5	*pscaling*ECON_SANC; % economic sanctions	
	
	% B_S
	u_bsbe = 4	*pscaling; % benefits of economic strength

	u_bst = 3	*pscaling; % effects of stability
	u_bsp = 3	*pscaling; % support of the people
	u_bsm = 2	*pscaling; % effect of the military 
	
	u_bscs = 3	*pscaling; % competition with G_S

	% P
	u_pgs = 2	*pscaling; % support for the G_S
	u_pge = 3	*pscaling; % economic benefits from G_E
	u_pbs = 2	*pscaling; % support for the B_S
	u_pbe = 3	*pscaling; % economic benefits from B_E
	u_pm = 5	*pscaling; % people dislike intrusions
	
	u_pt = 4	*pscaling; % people desire stability above all

	% M (note that most values are lower due to coupling of terms)
	
	u_mp = 3	*pscaling; % support or hindrance from the people
	u_mcse = 2	*pscaling; % responds to the threat of B_S with response proportional to G_E and G_S
	u_mb = 2	*pscaling; % military suffers from B_S, scaled by B_E

	u_mg = 2	*pscaling; % support from the G_S
	u_mtgs = 4	*pscaling; % if stable, decrease military proportional to support; if unstable, increase proportional to support
	u_mk = 5	*pscaling*MIL_SANC; % military sanctions
	
	% T
	u_tgs = 3	*pscaling; % state brings stability
	u_tbs = 3	*pscaling; % B_S brings instability
	
	u_tge = 2	*pscaling; % economic strength means stability, regardless of the source
	u_tbe = 2	*pscaling; 

	u_tg = 3	*pscaling; % military brings stability, but is only as effective as the state that wields it; would be 4 but have to adjust for coupling
	u_tp = 3	*pscaling; % whom do the people support?
	u_tb = 4	*pscaling; % B_S wields instability just as G_S wields the military; would be 4 but have to adjust for ln
	
	u_tce = 1	*pscaling; % negative effects from economic competition
	u_cex = 1	*pscaling; % adjust for the exponential
	
	u_tcs = 2	*pscaling; % negative effects from state competition
	u_csx = 1	*pscaling; % adjust for the exponential

for n=1:numiter
	G_E_n = G_E(n) * exp(u_gegs*G_S(n)  ...
		+ u_get*T(n) - u_gem*M(n)...
		- u_gece*B_E(n)) ...
		- u_gek*G_E(n);

	G_S_n = G_S(n) * exp(u_gsge*G_E(n) + u_gsm*M(n) ...
		+ u_gsp*P(n) + u_gst*T(n) ...
		- u_gscs*B_S(n));

	B_E_n = B_E(n) * exp(u_bebs*B_S(n) ...
		 - u_bet*T(n) ...
		 - u_becs*G_S(n) - u_bece*G_E(n)) ...
		 - u_bek*B_E(n);

	B_S_n = B_S(n) * exp(u_bsbe*B_E(n)  ...
		- u_bst*T(n) - u_bsp*P(n) - u_bsm*M(n)...
		- u_bscs*G_S(n));

	P_n = P(n) + u_pt*T(n) ...
		+ u_pgs*G_S(n) + u_pge*G_E(n) ...
		- u_pbs*B_S(n) - u_pbe*B_E(n)...
		- u_pm * M(n);
	
	M_n = M(n) * exp(-u_mg*G_E(n)*log(G_S(n)) + u_mcse*B_S(n)*G_S(n)*G_E(n)...
		+ u_mp*P(n) - u_mb*B_S(n)*B_E(n)...
		- u_mtgs*T(n)*G_S(n)*G_E(n))...
		- u_mk*M(n);

	T_n = T(n) + u_tgs*G_S(n) - u_tbs*B_S(n) + u_tge*G_E(n) + u_tbe*B_E(n)...
		+ u_tp*P(n) + u_tg*G_S(n)*M(n) + u_tb*B_E(n)*log(B_S(n)) ...
		- u_tcs*exp(u_csx*G_S(n)*B_S(n)) - u_tce*exp(u_cex*G_E(n)*B_E(n));

	if ENABLE_BRAKES
		if G_S_n > HUGE_VAL || B_S_n > HUGE_VAL
			disp(sprintf('BRAKED at %d',n));
			break;
		end
	end
	
	G_E(n+1) = G_E_n;
	G_S(n+1) = G_S_n;
	B_E(n+1) = B_E_n;
	B_S(n+1) = B_S_n;
	P(n+1) = P_n;
	M(n+1) = M_n;
	T(n+1) = T_n;
end

pstart = 1;
pend = n;

if PLOT_SUBSET
	pend = min(pss_end,n);
	pstart = max(pss_start,1);
end

figure;
plot(pstart:pend,G_E(pstart:pend),'c-');
hold on;
plot(pstart:pend,G_S(pstart:pend),'b-');
hold on;
plot(pstart:pend,B_E(pstart:pend),'y-');
hold on;
plot(pstart:pend,B_S(pstart:pend),'r-');
hold on;
plot(pstart:pend,P(pstart:pend),'m-');
hold on;
plot(pstart:pend,M(pstart:pend),'g-');
hold on;
plot(pstart:pend,T(pstart:pend),'w-');
legend('G_E','G_S','B_E','B_S','P','M','T');
% set(gca,'XColor',[11/255 132/255 199/255]); % useful for saving plots
% set(gca,'YColor',[11/255 132/255 199/255]); % and still see the axes
xlabel('n');
ylabel('arbitrary units');