/*
Author: Klaus Reiner Schenk-Hoppé (klaus@iew.unizh.ch)
        http://www.iew.unizh.ch/home/klaus
Last update: February 1, 2000

gcc -lm dyn_ineff_cond.c
*/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

/***  global variables  ***/

FILE            *fp_mat;


double          preiterations, iterations, steps;
double          stoch, fcase;
double          smin,smax;
double          alpha, A;
double          xi,xi1,xi2;
double          delta, delta1, delta2;
double          n, n1, n2;
double          *s, *Et;
/* double          test;    TMP */


double k_cd(double k, int i) {
       return ((1.0-delta) * k + s[i] * xi * A * pow(k,alpha))/(1.0+n);
}

double cond_cd(double k) {
      return (delta + n - xi * A * alpha * pow(k,alpha-1));
}

double k_ces(double k, int i) {
      return ((1.0-delta)*k+s[i]*xi*pow(1-A+A*pow(k,alpha),1.0/alpha))/(1.0+n);
}

double cond_ces(double k) {
      return  (delta + n - xi * A * pow(k,alpha-1.0) * pow(1-A+A*pow(k,alpha),1.0/alpha-1.0));
}


void Econd_cd(int i) {

  unsigned long  j;
  double k,c;
  double z;

  k = 100.0;

  Et[i] = 0.0;

  if (stoch == 0){
    xi = 1;
    n = 0;
    delta = 0.5;
    for(j=1;j<preiterations;j++)
     k = k_cd(k,i);
    for(j=1;j<iterations;j++) {
      Et[i] = (cond_cd(k) + (j-1)*Et[i])/j;
      k = k_cd(k,i);
    }
 }
 else{
   if (stoch == 1){
    n = 0;
    delta = 0.5;
   for(j=1;j<preiterations;j++) {
     z=drand48();
     if (z > 0.9){
         if (xi != xi1)
           xi = xi1;
       else
           xi = xi2;
     }
    k = k_cd(k,i);
  }
  for(j=1;j<iterations;j++) {
     Et[i] = (cond_cd(k) + (j-1)*Et[i])/j;
     z=drand48();
     if (z > 0.9){
         if (xi != xi1)
           xi = xi1;
       else
           xi = xi2;
     }
    k = k_cd(k,i);
  }
  }

  if (stoch == 2){
    xi = 1;
    delta = 0.5;
   for(j=1;j<preiterations;j++) {
     z=drand48();
     if (z > 0.9){
       if (n != n1)
            n = n1;
        else
            n = n2;
     }
    k = k_cd(k,i);
  }
  for(j=1;j<iterations;j++) {
    Et[i] = (cond_cd(k) + (j-1)*Et[i])/j;
    z=drand48();
    if (z > 0.9){
      if (n != n1)
           n = n1;
       else
           n = n2;
    }
    k = k_cd(k,i);
  }
  }
  if (stoch == 3){
    xi = 1;
    n = 0;
   for(j=1;j<preiterations;j++) {
     z=drand48();
     if (z > 0.9){
        if (delta != delta1)
            delta = delta1;
        else
            delta = delta2;
     }
    k = k_cd(k,i);
  }
  for(j=1;j<iterations;j++) {
    Et[i] = (cond_cd(k) + (j-1)*Et[i])/j;
    z=drand48();
    if (z > 0.9){
        if (delta != delta1)
            delta = delta1;
        else
            delta = delta2;
    }
    k = k_cd(k,i);
  }
  }
  }
}


void Econd_ces(int i) {

  unsigned long  j;
  double k,c;
  double z;

  k = 2.0;

  Et[i] = 0.0;

  if (stoch == 0){
    xi = 1;
    n = 0;
    delta = 0.5;
  for(j=1;j<preiterations;j++)
    k = k_ces(k,i);
  for(j=1;j<iterations;j++) {
    Et[i] = (cond_ces(k) + (j-1)*Et[i])/j;
    k = k_ces(k,i);
  }
 }
 else{
   if (stoch == 1){
    n = 0;
    delta = 0.5;
   for(j=1;j<preiterations;j++) {
     z=drand48();
     if (z > 0.9){
         if (xi != xi1)
           xi = xi1;
       else
           xi = xi2;
     }
    k = k_ces(k,i);
  }
  for(j=1;j<iterations;j++) {
    Et[i] = (cond_ces(k) + (j-1)*Et[i])/j;
     z=drand48();
     if (z > 0.9){
         if (xi != xi1)
           xi = xi1;
       else
           xi = xi2;
     }
    k = k_ces(k,i);
  }
  }

  if (stoch == 2){
    xi = 1;
    delta = 0.5;
   for(j=1;j<preiterations;j++) {
     z=drand48();
     if (z > 0.9){
       if (n != n1)
            n = n1;
        else
            n = n2;
     }
    k = k_ces(k,i);
  }
  for(j=1;j<iterations;j++) {
    Et[i] = (cond_ces(k) + (j-1)*Et[i])/j;
    z=drand48();
    if (z > 0.9){
      if (n != n1)
           n = n1;
       else
           n = n2;
    }
    k = k_ces(k,i);
  }
  }
  if (stoch == 3){
    xi = 1;
    n = 0;
   for(j=1;j<preiterations;j++) {
     z=drand48();
     if (z > 0.9){
        if (delta != delta1)
            delta = delta1;
        else
            delta = delta2;
     }
    k = k_ces(k,i);
  }
  for(j=1;j<iterations;j++) {
    Et[i] = (cond_ces(k) + (j-1)*Et[i])/j;
    z=drand48();
    if (z > 0.9){
        if (delta != delta1)
            delta = delta1;
        else
            delta = delta2;
    }
    k = k_ces(k,i);
  }
  }
  }
}




typedef struct {
     long type;   /* type */
     long mrows;  /* row dimension */
     long ncols;  /* column dimension */
     long imagf;  /* flag indicating imag part */
     long namlen; /* name length (including NULL) */
} Fmatrix;



/*
 * savemat - C language routine to save a matrix in a MAT-file.
 *
 * Here is an example that uses 'savemat' to save two matrices to disk,
 * the second of which is complex:
 *
 *  FILE *fp;
 *  double xyz[1000], ar[1000], ai[1000];
 *  fp = fopen("foo.mat","wb");
 *  savemat(fp, 2000, "xyz", 2, 3, 0, xyz, (double *)0);
 *  savemat(fp, 2000, "a", 5, 5, 1, ar, ai);
 *      fclose(fp);
 *
 * Author J.N. Little 11-3-86
 */

void savemat(FILE *fp, int type, char *pname, int mrows, int ncols,
         int imagf, double *preal, double *pimag)
     /* FILE *fp;       File pointer */
     /* int type;       Type flag: Normally 0 for PC, 1000 for Sun, Mac, and     */
     /*                 Apollo, 2000 for VAX D-float, 3000 for VAX G-float    */
     /*                 Add 1 for text variables.    */
     /*                 See LOAD in reference section of guide for more info. */
     /* char *pname;    pointer to matrix name */
     /* int mrows;      row dimension */
     /* int ncols;      column dimension */
     /* int imagf;  imaginary flag */
     /* double *preal;  pointer to real data */
     /* double *pimag;  pointer to imag data */
{
    Fmatrix x;
    int mn;

    x.type = type;
    x.mrows = mrows;
    x.ncols = ncols;
    x.imagf = imagf;
    x.namlen = strlen(pname) + 1;
    mn = x.mrows * x.ncols;

    fwrite(&x, sizeof(Fmatrix), 1, fp);
    fwrite(pname, sizeof(char), x.namlen, fp);
    fwrite(preal, sizeof(double), mn, fp);
    if (imagf) {
         fwrite(pimag, sizeof(double), mn, fp);
    }
}



int main(int argc, char *argv[]) {

  unsigned long   i;
  double ds;

if (argc < 8) {
  fprintf(stderr,  "parameters: production-function  s_min  s_max  steps  stochastic-case  pre-iterations  iterations\n
                    production-function (1=Cobb-Douglas, 2=CES)\n
                    s_min  s_max  steps  (saving rate takes values in [s_min, s_max] on an equidistant grid with steps+1 vertices)\n
                    stochastic-case (0=deterministic, 1=production shocks (xi), 2=stoch. population growth rate (n), 3=stoch. depreciation (delta))\n
                    pre-iterations (recommended > 100)\n
                    iterations (recommended > 10,000) \n");
                    fflush(stderr);
                    exit(1);
}
else{
  sscanf(argv[1],  "%lf",  &fcase);
  sscanf(argv[2],  "%lf",  &smin);
  sscanf(argv[3],  "%lf",  &smax);
  sscanf(argv[4],  "%lf",  &steps);
  sscanf(argv[5],  "%lf",  &stoch);
  sscanf(argv[6],  "%lf",  &preiterations);
  sscanf(argv[7],  "%lf",  &iterations);
}


  s = malloc(sizeof(double)*(steps+1));
  Et = malloc(sizeof(double)*(steps+1));

  ds = (smax-smin)/((double)steps);

  if (fcase<2) {
    A = 10.0;
    alpha = 0.75;}
  else {
    A = 0.5;
    alpha = 0.5;
  }

  xi1 = 0.75;
  xi2 = 1.25;
  xi = xi2;

  n1 = -0.07;
  n2 =  0.07;
  n = n2;

  delta1 = 1.0/3.0;
  delta2 = 2.0/3.0;
  delta = delta2;


if (fcase<2) {
  for (i=0; i < steps+1; i++) {
   s[i] = smin + i * ds;
   Econd_cd(i);
   fprintf(stderr,"s = %f \t Et = %f \n",s[i], Et[i]);
  }
}
else{
  for (i=0; i < steps+1; i++) {
   s[i] = smin + i * ds;
   Econd_ces(i);
   fprintf(stderr,"s = %f \t Et = %f \n",s[i], Et[i]);
  }
}




  if ((fp_mat = fopen("dyn_ineff_cond.mat","wb")) == NULL) {
    fprintf(stderr,"Unable to open output file dyn_ineff_cond.mat\n");
    fflush(stderr);
    exit(1);
  }


 savemat(fp_mat, 1000, "production-function", 1, 1, 0, &fcase, (double *)0);
 savemat(fp_mat, 1000, "alpha", 1, 1, 0, &alpha, (double *)0);
 savemat(fp_mat, 1000, "A", 1, 1, 0, &A, (double *)0);
 savemat(fp_mat, 1000, "preiterations", 1, 1, 0, &preiterations, (double *)0);
 savemat(fp_mat, 1000, "iterations", 1, 1, 0, &iterations, (double *)0);
 savemat(fp_mat, 1000, "s", 1, steps+1, 0, s, (double *)0);
 savemat(fp_mat, 1000, "Et", 1, steps+1, 0, Et, (double *)0);

 fclose(fp_mat);

  return 0;
}

