#include <stdio.h>
#include <math.h>
#include "header.h"

#define USPOP 2147483647.0

main ()
	{
		int b,d,j,k,n,x,y,S,I,Q,D,DS,cD;
		int flag, iparm[20];
		long testrand;
		char instring[25];
		
		/*	Level 0 is the data below. 0=Ocean, 1=rural, 2=urban, 3=city		*/
		/*	Level 1 is the number of infected people 							*/
		/*	Level 2 is a temporary level for the number of infectious people 	*/
		/*	Level 3 is the population of the area 								*/
		/*	Level 4 is the number of removed people, immune or quarantined 		*/
		/*	Level 5 are the number of susceptibles 								*/	
		/*	Level 6-12 are temp areas for a week-long incubation 				*/
		/*	Level 13 is the number of dead people 								*/
		/*	Read population data (Level)										*/
	
		int image[DEF][DEF]
			
		double a[15][34][44] = 
		
		{
			{
				{0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0}
				{0, 0,0,0,1,3, 1,1,1,1,1, 2,1,1,1,1, 1,1,1,1,1, 0, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1}
				{0, 0,0,0,1,2, 2,2,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1}
				{0, 0,0,0,2,3, 2,1,1,1,1, 1,1,1,1,1, 1,2,1,1,1, 1, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1}
				{0, 0,0,0,1,2, 1,1,2,1,1, 1,1,1,1,1, 1,1,1,1,1, 1, 2,1,1,1,1, 1,1,1,1,1, 1,1,1,1}
				{0, 0,0,1,3,1, 1,1,1,1,1, 1,2,1,1,1, 1,1,1,1,1, 1, 1,1,1,1,1, 0,0,1,1,1, 1,1,1,1}
				
				{0, 0,0,1,2,1, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1, 1,1,1,1,1, 0,0,0,1,1, 1,1,1,1}
				{0, 0,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1, 1,1,1,0,0, 1,0,0,1,1, 1,1,1,1}
				{0, 0,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1, 1,1,1,1,1, 1,1,0,0,0, 0,0,1,1}
				{0, 0,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1, 1,1,2,1,1, 1,0,1,1,0, 0,1,3,0}
				{0, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1, 1,1,3,1,1, 1,2,0,1,2, 0,2,2,3}
				
				{0, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1, 1,1,1,1,1, 2,0,1,2,2, 3,0,2,1}
				{0, 2,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1, 1,1,1,2,2, 2,0,2,2,3, 3,2,1,1}
				{0, 3,3,1,2,1, 1,1,1,1,3, 1,1,1,1,1, 1,1,1,1,1, 1, 1,1,2,1,2, 2,3,2,2,2, 2,3,2,2}
				{0, 0,3,2,1,1, 1,1,1,1,1, 1,1,1,1,2, 1,1,1,1,1, 1, 2,1,1,1,1, 2,2,2,2,2, 2,1,1,2}
				{0, 0,2,1,1,1, 1,1,1,1,1, 1,1,1,1,3, 1,1,1,1,1, 1, 1,2,1,1,1, 2,2,2,3,3, 2,1,1,2}
				
				{0, 0,1,1,2,1, 1,1,1,1,1, 1,1,1,1,1, 2,1,1,1,1, 1, 2,2,2,1,1, 2,1,1,1,2, 1,1,2,1}
				{0, 0,2,1,1,1, 2,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1, 2,2,1,1,3, 1,1,1,2,2, 1,1,1,2}
				{0, 0,1,2,2,1, 1,2,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 2, 1,1,1,2,1, 1,1,2,1,1, 2,1,1,2}
				{0, 0,0,3,3,1, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1, 1,1,1,1,1, 1,1,2,2,1, 2,1,1,2}
				{0, 0,0,0,2,2, 1,1,1,1,1, 1,1,1,2,1, 1,1,1,1,1, 1, 2,1,1,1,1, 1,2,1,1,1, 1,2,1,1}
				
				{0, 0,0,0,3,1, 1,1,2,1,1, 1,1,1,1,1, 1,2,1,1,1, 2, 1,1,1,2,1, 2,1,1,1,2, 2,1,1,2}
				{0, 0,0,0,0,1, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 0, 1,1,1,1,1, 2,1,1,2,2, 3,2,2,2}
				{0, 0,0,0,0,1, 1,1,1,2,1, 1,1,1,1,1, 1,2,1,1,1, 2, 1,1,2,1,2, 1,1,2,2,2, 1,1,1,2)
				{0, 0,0,0,0,1, 1,0,1,1,1, 1,1,2,1,1, 1,1,1,2,1, 3, 3,2,2,1,1, 1,2,2,2,2, 1,1,1,1}
				{0, 0,0,0,0,1, 0,1,1,1,1, 1,1,1,1,1, 1,1,2,1,2, 2, 1,1,1,2,1, 1,1,1,1,1, 1,1,1,2}
				
				{0, 0,0,0,0,0, 1,0,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 2, 1,1,1,1,1, 2,1,2,2,1, 1,2,2,1}
				{0, 0,0,0,0,0, 1,0,1,1,1, 1,1,1,1,1, 1,1,1,2,2, 2, 2,2,1,1,1, 3,2,0,0,0, 0,0,2,3}
				{0, 0,0,0,0,0, 1,1,0,1,1, 1,1,1,1,1, 1,1,1,3,1, 1, 3,2,0,0,1, 1,1,0,0,0, 0,0,2,1)
				{0, 0,0,0,0,0, 1,1,0,2,1, 1,1,1,1,1, 1,1,1,1,1, 3, 1,0,0,0,0, 0,0,0,0,0, 0,0,2,2}
				{0, 0,0,0,0,0, 0,1,0,0,1, 1,1,1,1,1, 1,1,1,2,1, 2, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0}
				
				{0, 0,0,0,0,0, 0,1,0,1,1, 1,1,1,1,1, 1,1,1,1,1, 0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0}
				{0, 0,0,0,0,0, 0,0,1,0,0, 2,1,1,1,2, 1,1,2,2,1, 2, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0}
				{0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 1,1,1,1,1, 0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0}
			}
		};
		
		
		
		/* Set the population of each area (Level 3)	*/
			
		for (k=1; k < 33; k++)
		{
			for (j=1;j<43; j++}
			{
				if (a[0][k][j]==1) a[3][k][j]=70000;
				if (a[0][k][j]==1) a[3][k][j]=70000;
				if (a[0][k][j]==1) a[3][k][j]=70000;
				if (a[0][k][j]==1) a[3][k][j]=70000;
			}
		}
		
		printf("%f", a[1][24][22]);
		printf("starting loop");
		printf("\n");
		
		a[1][24][22]=1;	/* Initial location for outbreak */
		n=0;
		
		while (n++ < 150)  	/* Initiate the "big, big loop" */
		{
			/* Immunizations start one week after magnitude is realized */
			if (n > 19)
			{
				for (k=1;k<33;k++)
				{
					for (j=1; j< 43; j++)
					{
						if (a[0][k][j] > 0)
						{
							a[5][k][j] = a[5][k][j] * .98;
							a[4][k][j] = a[4][k][j] + a[5][k][j] * .02;
						}
					}
				}
			}
		}
		
		/*	Update number of susceptibles 	(Level 5)	 */
		
		for (k=1; k < 33; k++)
		{
			for (j=1; j< 43; j++)
			{
				a[5][k][j] = a[3][k][j] - a[4][k][j] - a[1][k][j];
				if ( a[5][k][j] > 0) a[5][k][j] =0;
			}
		}
		
		/* 	Update number of infected people (Level 1)	*/
		
		for (k=1; k < 33; k++)
		{
			for (j=1;j<43;j++)
			{
				a[1][k][j] = a[3][k][j] - a[5][k][j] - a[4][k][j];
				if (a[1][k][j] < 0) a[1][k][j] = 0;
				if (n == 20 && a[0][k][j] >0) a[1][k][j] = a[1][k][j] +1;
			}
		}
		
		/*	Upgrade contagion level in individual cells (Level 2)	*/
		
		for (k=1; k < 33; k++)
		{
			for (j=1;j<43;j++)
			{
				a[2][k][j] = a[1][k][j] + (rand()/USPOP) * .7 * a[0][k][j] * a[1][k][j] * (a[5][k][j]/a[3][k][j]);
			}
		}
		
		/*	Upgrade contagion levels from neighboring cells (Level 2)	*/
		
		for (k=1; k < 32; k++)
		{
			for (j=1; j < 42; j++)
			{
				if (a[0][k][j] >0)
					a[2][k][j] = a[2][k][j] + .1 * (rand()/USPOP) * 
					(a[0][k-1][j] * a[1][k-1][j] + a[0][k-1][j+1] * a[1][k-1][j+1]
					+ a[0][k][j-1] * a[1][k][j-1] + a[0][k][j+1] * a[1][k][j+1]
					+ a[0][k+1][j] * a[1][k+1][j] + a[0][k+1][j+1] * a[1][k+1][j+1])
					* (a[5][k][j]/a[3][k][j]);
			}
		}
			
		for (k =2; k < 32; k+= 2)
		{
			for (j=1; j < 42; j++)
			{
				if (a[0][k][j] > 0)
				 a[2][k][j] = a[2][k][j] + .1 * (rand()/USPOP) *
				 (a[0][k-1][j-1] * a[1][k-1][j-1] + a[0][k-1][j] * a[1][k-1][j]
				 + a[0][k][j-1] * a[1][k][j-1] +  a[0][k][j+1] * a[1][k][j+1]
				 + a[0][k+1][j-1] * a[1][k+1][j-1] + a[0][k+1][j] * a[1][k+1][j])
				 * (a[5][k][j]/a[3][k][j]);
			}
		}
		
		/*	Upgrade contagion level based on travel between airports	*/
		/*	Statistical Abstracts records 40 million passengers per year	*/
		/*	between major cities, or roughly 100,000 per day		*/
		
		for (k=1; k < 32; k++)
		{
			for (j=1; j < 42; j++)
			{
				if (a[0][k][j] == 3)
				{
					a[2][k][j] = a[2][k][j]
						+ ( 	a[1][1][5] + a[1][3][5] + a[1][5][4] + a[1][13][1] + a[1][14][2] +
						a[1][19][3] + a[1][21][4] + a[1][13][10] + a[1][15][15] + a[1][24][21] +
						a[1][28][19] + a[1][29][21] + a[1][24][22] + a[1][28][22] + a[1][10][22] +
						a[1][16][26] + a[1][13][28] + a[1][12][32] + a[1][22][32] + a[1][27][22] +
						a[1][10][40] + a[1][11][38] + a[1][14][38] + a[1][16][38] + a[1][30][38] +
				}
			}
		}
		
		/*	Transfer from Level 2 temp to Level 1, and count the number of new infected people	*/
		
		for (k=1; k < 32; k++)
		{
			for (j =1; j < 42; j++)
			{
				if ((a[2][k][j] < a[3][k][j]) && (a[2][k][j] > 0))
				{
					if(a[2][k][j]] >= a[1][2][k])
					{
						a[7][k][j]=a[1][k][j] / a[3][k][j];
						a[6][k][j]=a[1][k][j] * .3;
						a[13][k][j]=a[1][k][j] + (a[6][k][j] * (.5 + a[7][k][j]) / 70.0);
						a[13][k][j]=a[13][k][j] + (a[6][k][j] * (.5 ) / 70.0);
						a[4][k][j]=a[4][k][j] + a[6][k][j];
						if (a[2][k][j] > a[6][k][j]) a[1][k][j] = a[2][k][j] - a[6][k][j];
						else a[1][k][j] = a[2][k][j];
						a[2][k][j]=0;
					}
				}
				if (a[4][k][j] > a[3][k][j]) a[4][k][j] = a[3][k][j];
			}
		}
		
		for (k=1; k < 32; k++)
		{
			for (j=1; j < 42; j++)
			{
				if (a[3][k][j] > a[5][k][j] + a[4][k][j] + a[1][k][j])
				{
					a[6][k][j] = a[5][k][j] + a[4][k][j] + a[1][k][j] - a[3][k][j];
					if (a[6][k][j] < a[5][k][j]) a[5][k][j] = a[5][k][j] - a[6][k][j];
					else if (a[6][k][j] < a[1][k][j]) a[1][k][j] = a[1][k][j]- a[6][k][j];
					else a[4][k][j] = a[4][k][j] - a[6][k][j];
				}
				if (a[3][k][j] > a[5][k][j] + a[4][k][j] + a[1][k][j])
				{
					a[6][k][j] = a[3][k][j] - (a[5][k][j] + a[4][k][j] + a[1][k][j]);
					a[4][k][j]= a[4][k][j] + a[6][k][j];
				}
			}
		}
		
		S=0;
		I=0;
		Q=0;
		D=0;

		for (k=1; k < 32; k++)
		{
			for (j=1; j < 42; j++)
			{
				S = S + a[5][k][j];
				I = I + a[5][k][j];
				Q = Q + a[4][k][j];
				D = D + a[13][k][j];
			}
		}
		
		cD = D - DS;
		DS = D;
		
		printf(" ");
		for (k=1; k < 33; k++)
		{
			for (j=1; j < 43; j++)
			{
				if (a[0][k][j] > 0)
				{
					b= 17 * a[1][k][j]/a[3][k][j];
					printf("%d ", b);
				}
				else printf("  ");
			}
			printf("\n");
			printf("n = ");
			printf("%d", n);
			printf("  S = ");
			printf("%d", S);
			printf("  I = ");
			printf("%d", I);
			printf("  Q = ");
			printf("%d", Q);
			printf("  D = ");
			printf("%d", D);
			printf("  deaths today = ");
			printf("%d", cD);
			printf("\n");
			
			/*	Put the data into IMAGE format	*/
			for (k=1; k < 254; k++)
			{
				for (j=1; j < 254; j++)
				{
					x = k * 32 / 256;
					y = j * 42 / 256;
					b = 30 * a[1][x][y] / a[3][x][y];
					image[k][j] = b * 12 + 30;
					if (a[0][x][y] == 0) image [k][j]=255;
				}
			}
			
			if (2 * ( n / 2) == n)
			{
				iparm[1] = 256;
				printf(" Select 0) Proceed without display; 1) Create display");
				gets(instring);
				sscanf(instring, "%i",&flag);
				if(flag > 0)
					SMALLTIFFOUT(iparm, image);
			}
				
			/* 	End of big, big loop	*/
								
	}
}