#include <stdio.h>
#include <iostream.h>
#include <fstream.h>
#include <math.h>

//get slope function
float get_slope(float y, float x, float c)

{
float m;
m = (y-c)/(x-1);
return(m);
}



//main program
int main()
{

//declare internal variables
int age, d, choice;


//declare life-history variables
int L2span, adspan, respan, glife, gopen;
float adjbrate, adults, aphids, soldiers, youngvir, youngre;
float E, k, a, Pmax, soldmort, Esold, Eother;
float brate, brateini, bratefin, brateincr;
float sfract, sfractini, sfractfin, sfractincr;
float qini, qfin, qincr;
int q, qlong;


//open input and output files
char fileout[35], filein[35];
cout << "Summary output (1) or day-by-day tracking (2)?";
cin >> choice;
cout << "Please select input filename\n";
cin >> filein;
cout << "Please select output filename\n";
cin >> fileout;

ofstream fout(fileout);
ifstream fin(filein);
char params[201];
fin.getline (params, 200);

if (choice == 1) {
   fout << params<<"\t"<<"incr aphids\t"<<"peak aphids\t"<<"peak soldiers\t"<<"peak youngvir\t"<<"peak virgins\t"<<"peak youngre\t"<<"peak reproductives\t";
   fout << "\t"<<"aphids\t"<<"soldiers\t"<<"youngvir\t"<<"virgins\t"<<"youngre\t"<<"reproductives\n";
   }

if (choice == 2) {
   fout << params<<"\n";
   }


//cycle until no more data in the file
glife = 1;
while (glife >0) {
fin >> glife>>qini>>qfin>>L2span>>adspan>>respan>>brateini>>bratefin>>sfractini>>sfractfin>>Pmax>>k>>a>>gopen>>soldmort;

if (glife>0) {
   fout << glife<<"\t"<<qini<<"\t"<<qfin<<"\t"<<L2span<<"\t"<<adspan<<"\t"<<respan<<"\t"<<brateini<<"\t"<<bratefin<<"\t"<<sfractini<<"\t"<<sfractfin<<"\t"<<Pmax<<"\t"<<k<<"\t"<<a<<"\t"<<gopen<<"\t"<<soldmort;
   if (choice == 2) {
      fout << "\n";
      fout << "day"<<"\t"<<"aphids"<<"\t"<<"soldiers"<<"\t"<<"youngvir\t"<<"virgins\t"<<"youngre\t"<<"reproductives\n";
      }
   }


//calculate longest L1 over season
if  (qini > qfin)
    qlong = int(qini);
else
    qlong = int(qfin);


//zero arrays and starting position
float L1[1000], L2[1000], ad[1000], re[1000],
 readults, die_of_old_age;
for (age=0; age<(qlong+1); age++) {
           L1[age] = 0;
           }
for (age=0; age<(L2span+1); age++) {
           L2[age] = 0;
           }
for (age=0; age<(adspan+1); age++) {
           ad[age] = 0;
           }
for (age=0; age<(respan+1); age++) {
           re[age] = 0;
           }
           
readults = 0;
ad [1] = 10;
Eother = Esold = 0;
qincr = get_slope(qfin, glife, qini);
sfractincr = get_slope(sfractfin, glife, sfractini);
brateincr = get_slope(bratefin, glife, brateini);


//iteration for each day of the gall
  for (d=1; d<(glife+1); d++) {


       //calculate number of adults giving birth and total number of aphids
       adults = 0;
       aphids = 0;
       soldiers = 0;
       youngvir = 0;
       youngre = 0;
       for (age=1; age<(qlong+1); age++) {
           aphids += L1[age];
           soldiers += L1[age];
           }
       for (age=1; age<(L2span+1); age++) {
           aphids += L2[age];
           youngvir += L2[age];
           }
       for (age=1; age<(adspan+1); age++) {
           aphids += ad[age];
           adults += ad[age];
           }
       for (age=1; age<(respan+1); age++) {
           aphids += re[age];
           youngre += re[age];
           }


       //calculate mortality/predation
       if  ((gopen==1)||((gopen==0)&&(readults>0))) {
           E = Pmax*((1-exp(-k*aphids))/(aphids + a * soldiers));
           if (soldiers>0)
              Esold = (E * aphids)*((soldmort*soldiers)/(aphids-soldiers+soldmort*soldiers))/soldiers;
           else
               Esold = 0;
           Eother = (E *aphids)*((aphids-soldiers)/(aphids-soldiers+soldmort*soldiers))/(aphids-soldiers);
           }
       else
           E = 0;
           
           
       //calculate soldier fraction and soldier instar duration
       sfract = (sfractincr*d+sfractini-sfractincr);
       q = int(qincr*d+qini-qincr);
       
       
       //calculate birthrate and its density-dependent adjustment
       brate = (brateincr*d+brateini-brateincr);
       if ((1-aphids/1000)>0)
          adjbrate = brate*(1-aphids/1000);
       else
          adjbrate = 0;
           
           
	   //old adults die
	   die_of_old_age = ad[adspan];
	   
       //calculate new numbers of adults
       for (age=(adspan-1); age>0; age--) {
           ad[age+1]=ad[age]*(1-Eother);
           }


       //old L2 move to adults
       ad[1]=L2[L2span]*(1-Eother);

       //calculate new numbers of L2
       for (age=(L2span-1); age>0; age--) {
           L2[age+1] = L2[age]*(1-Eother);
           }


       //old L1 move to L2
       L2[1] = L1[q]*(1-Eother);
       for (age=qlong; age>q; age--) {
           L2[1] += L1[age]*(1-Eother);
           L1[age] = 0;
           }

       //calculate new numbers of L1
       for (age=(q-1); age>0; age--) {
           L1[age+1] = L1[age]*(1-Esold);
           }


       //calculate number of immature reproductives maturing and leaving gall
       readults = readults+re[respan]*(1-Eother);

       //calculate new numbers of immature reproductives
       for (age=(respan-1); age>0; age--) {
           re[age+1] = re[age]*(1-Eother);
           }


       //new soldiers and immature reproductives are born
       L1[1] = adults*adjbrate*sfract*(1-Esold);
       re[1] = adults*adjbrate*(1-sfract)*(1-Eother);


	   //output
       if (choice == 1){
          if (d == 50)
             fout << "\t"<<aphids;
          if (d == 100)
              fout << "\t"<<aphids<<"\t"<<soldiers<<"\t"<<youngvir<<"\t"<<adults<<"\t"<<youngre<<"\t"<<readults<<"\t";
       }

       if (choice == 2)
          fout << d<<"\t"<<aphids<<"\t"<<soldiers<<"\t"<<youngvir<<"\t"<<adults<<"\t"<<youngre<<"\t"<<readults<<"\n";
       }
       
    if (glife>0){
       if (choice == 1)
       fout << "\t"<<aphids<<"\t"<<soldiers<<"\t"<<youngvir<<"\t"<<adults<<"\t"<<youngre<< "\t"<<readults<<"\n";
       }
  }


cout << "end\n";
return(0);
}