void main ()
{

/* information for model and integration of ODEs*/
/* new statements and braces*/
float rs=0.5;
while (ENDTIME<= 160000000)
{ rs=0.0;
}
/* end of new statements and braces, except rs=0.5 is removed from below float statement*/

float ct, st, s2t, 
dmp, dp0, dp1, dp2, dpn, mp2, p02, p12, p22, pn2, 
fi=0.0, tytime_on=0.0, timer_off=0.0, t_ran_on=0.0, t_ran_off=0.0, kran, 
r1=6.0, r2=3.0, r3=6.0, r4=3.0, rd=1.5, rm=0.3 ;

/* initial concetrations nM */

float mp, p0, p1, p2, pn, max1=1.0, max2=1.0, max3=1.0, mp1=2.0,
		pn1=2.2,
		p01=2.5,
		p11=2.1,
		p21=2.5,
		mpij, mp12, mp23, mp34, mp45,
		pnij, pn12, pn23, pn34, pn45,
		p0ij, p012, p023, p034, p045,
		p1ij, p112, p123, p134, p145,
		p2ij, p212, p223, p234, p245;

int t=0, sch=1, ITERATION=1;
float dt=1.0/BETA;			/* this makes dt = 0.00005*/

long id = PERT_SEED;
long *idum;
float sign=1.0;

char filename[300];
FILE *gg;

	gg = fopen("22mar07v2.dat", "w");		/*change filename here*/

idum=&id;

/* Time Loop */

t=0;

/*square wave generator*/

while (t++ < ENDTIME)
{
	if (tytime_on <= (12.0 + t_ran_on*bandwidth) && sch == 1)
		{
		fi=mp_max;
		tytime_on = tytime_on + dt;
		}
		else if (tytime_on > (12.0 + t_ran_on*bandwidth) && sch==1)
		{
		tytime_on = 0.0;
		sch=2;
		t_ran_off = gasdev(idum);
		while (t_ran_off <-2.0 || t_ran_off > 2.0) t_ran_off = gasdev(idum);
		}
		
		if (timer_off <= (12.0 + t_ran_off*bandwidth) && sch == 2)
		{
		fi=mp_min;
		timer_off = timer_off + dt;
		}
		else if (timer_off > (12.0 + t_ran_off*bandwidth) && sch==2)
		{
		timer_off = 0.0;
		sch=1;
		t_ran_on = gasdev(idum);
		while (t_ran_on <-2.0 || t_ran_on > 2.0) t_ran_on = gasdev(idum);
		/* fprintf(stderr, "%f %f \n", kran, t_ran_on);*/
		}

/*step for Runge Kutta Method*/

mp=mp1;
pn=pn1;
p0=p01;
p1=p11;
p2=p21;

/*ODEs for step one of Runge Kutta Method*/

mpij = (rs+fi)*(pow(Ki,n))/(pow(Ki,n) + pow(pn, n)) - rm*(mp/(Km + mp));
pnij = k1*p2 - k2*pn;
p0ij = ks*mp - r1*(p0/(K1+p0)) + r2*(p1/(K2+p1));
p1ij = r1*(p0/(K1+p0)) - r2*(p1/(K2+p1)) - r3*(p1/(K3+p1)) + r4*(p2/(K4+p2));
p2ij = r3*(p1/(K3+p1)) - r4*(p2/(K4+p2)) - rd*(p2/(Kd+p2)) - k1*p2 + k2*pn;

/* printf ( "%f %f %f %f %f %f\n", t*dt, mpij*1000, pnij*1000, p0ij*1000, p1ij*1000, p2ij*1000); */

/*step for Runge Kutta Method*/

mp12=dt*mpij;
pn12=dt*pnij;
p012=dt*p0ij;
p112=dt*p1ij;
p212=dt*p2ij;

mp=mp1+mp12;
pn=pn1+pn12;
p0=p01+p012;
p1=p11+p112;
p2=p21+p212;

/*ODEs for step two of Runge Kutta Method*/

mpij = (rs+fi)*(pow(Ki,n))/(pow(Ki,n) + pow(pn, n)) - rm*(mp/(Km + mp));
pnij = k1*p2 - k2*pn;
p0ij = ks*mp - r1*(p0/(K1+p0)) + r2*(p1/(K2+p1));
p1ij = r1*(p0/(K1+p0)) - r2*(p1/(K2+p1)) - r3*(p1/(K3+p1)) + r4*(p2/(K4+p2));
p2ij = r3*(p1/(K3+p1)) - r4*(p2/(K4+p2)) - rd*(p2/(Kd+p2)) - k1*p2 + k2*pn;

/* printf ( "%f %f %f %f %f %f\n", t*dt, mpij*1000, pnij*1000, p0ij*1000, p1ij*1000, p2ij*1000); */

/*step for Runge Kutta Method*/

mp23=dt*mpij;
pn23=dt*pnij;
p023=dt*p0ij;
p123=dt*p1ij;
p223=dt*p2ij;

mp=mp1+mp23;
pn=pn1+pn23;
p0=p01+p023;
p1=p11+p123;
p2=p21+p223;

/*ODEs for step three of Runge Kutta Method*/

mpij = (rs+fi)*(pow(Ki,n))/(pow(Ki,n) + pow(pn, n)) - rm*(mp/(Km + mp));
pnij = k1*p2 - k2*pn;
p0ij = ks*mp - r1*(p0/(K1+p0)) + r2*(p1/(K2+p1));
p1ij = r1*(p0/(K1+p0)) - r2*(p1/(K2+p1)) - r3*(p1/(K3+p1)) + r4*(p2/(K4+p2));
p2ij = r3*(p1/(K3+p1)) - r4*(p2/(K4+p2)) - rd*(p2/(Kd+p2)) - k1*p2 + k2*pn;

/* printf ( "%f %f %f %f %f %f\n", t*dt, mpij*1000, pnij*1000, p0ij*1000, p1ij*1000, p2ij*1000); */

/*step for Runge Kutta Method*/

mp34=dt*mpij;
pn34=dt*pnij;
p034=dt*p0ij;
p134=dt*p1ij;
p234=dt*p2ij;

mp=mp1+mp34;
pn=pn1+pn34;
p0=p01+p034;
p1=p11+p134;
p2=p21+p234;

/*ODEs for step four of Runge Kutta Method*/

mpij = (rs+fi)*(pow(Ki,n))/(pow(Ki,n) + pow(pn, n)) - rm*(mp/(Km + mp));
pnij = k1*p2 - k2*pn;
p0ij = ks*mp - r1*(p0/(K1+p0)) + r2*(p1/(K2+p1));
p1ij = r1*(p0/(K1+p0)) - r2*(p1/(K2+p1)) - r3*(p1/(K3+p1)) + r4*(p2/(K4+p2));
p2ij = r3*(p1/(K3+p1)) - r4*(p2/(K4+p2)) - rd*(p2/(Kd+p2)) - k1*p2 + k2*pn;

/* printf ( "%f %f %f %f %f %f\n", t*dt, mpij*1000, pnij*1000, p0ij*1000, p1ij*1000, p2ij*1000);*/

/*step for Runge Kutta Method*/

mp45=dt*mpij;
pn45=dt*pnij;
p045=dt*p0ij;
p145=dt*p1ij;
p245=dt*p2ij;

/* Runge Kutta Method */

mp2=mp1 + (mp12 +2.0*mp23 +2.0*mp34 +mp45)/6.0 ;
pn2=pn1 + (pn12 +2.0*pn23 +2.0*pn34 +pn45)/6.0 ;
p02=p01 + (p012 +2.0*p023 +2.0*p034 +p045)/6.0 ;
p12=p11 + (p112 +2.0*p123 +2.0*p134 +p145)/6.0 ;
p22=p21 + (p212 +2.0*p223 +2.0*p234 +p245)/6.0 ;

if (mp2>LOW_LIMIT) mp1=mp2; else mp1=LOW_LIMIT;
if (pn2>LOW_LIMIT) pn1=pn2; else pn1=LOW_LIMIT;
if (p02>LOW_LIMIT) p01=p02; else p01=LOW_LIMIT;
if (p12>LOW_LIMIT) p11=p12; else p11=LOW_LIMIT;
if (p22>LOW_LIMIT) p21=mp2; else p21=LOW_LIMIT;

/* 	if(t>10000000) {
	if(!(t%2)) {
	max3=max2;
	max2=max1;
	max1=pn1;
	
	if (max2>max1 && max2>max3) {
	fprintf (gg, "%f %f\n, t_on, max2);
	} } } */
	
if(!(t%200000))
printf ("%f %f %f %f\n", t*dt, mp1, pn1, p01);

if (!(t%1000))
{
	fprintf (gg, "%f %f %f %f %f %f %f\n", t*dt, mp1, pn1, p01, p11, p21, fi);

}

}	/* end t */

	/* end main*/
}