#include <stdio.h>
#include <math.h>
#include "geostat.h"
#include "genlib.h"


/*addstat                                                 */
/*adds mean and variance effects to the N(0,1) realization*/
/*input:                                                  */
/*realin: structure defining a realization -              */
/*        must have zeio mean and unit variance           */
/*stat: structure defining the mean and variance          */
/*output:                                                 */
/*realout: structure defining a realization -             */
void addstat(struct realization_mod *realin,struct statistic_mod stat ,struct realization_mod *realout)

{

  int i,nblockm,nblockv;
  double r,moy,var;

  /*Is the output realization allocated ?*/
  /*is its length equal to realin.n?*/
  if ((*realout).vector == NULL || (*realout).n != (*realin).n) {
    (*realout).vector = (double *) malloc((*realin).n * sizeof(double));
    if ((*realout).vector == NULL)
      Error("No memory available");
  }
  (*realout).n = (*realin).n;


  /*test over the input realization*/
  switch ((*realin).code) {
  case 0:
  case 1:
    (*realout).code = 2;
    break;
  case 6:
    (*realout).code = 7;
    break;
  default:
    (*realout).code = (*realin).code;
    break;
  }

  for (i = 0; i < (*realin).n; i++) {


    /*mean*/
    switch (stat.nblock_mean) {
    case 1:
      /*stationary case*/
      nblockm = 1;
      break;
    default:
      /*number of the cell - nonstationary case*/
      nblockm = i+1;
      break;
    }


    /*variance*/
    switch (stat.nblock_var) {
    case 1:
      /*stationary case*/
      nblockv = 1;
      break;
    default:
      /*number of the cell - nonstationary case*/
      nblockv = i+1;
      break;
    }


    switch (stat.type) {
    case 0:
      /*normal case*/
      moy = stat.mean[nblockm-1];
      var = stat.variance[nblockv-1];
      break;
    case 1:
      /*lognormal (natural) case*/
      r = (double)sqrt(stat.variance[nblockv-1])/stat.mean[nblockm-1];
      r *= r;
      moy = (double)log(stat.mean[nblockm-1]/sqrt(1.0+r));
      var = (double)log(1.0+r);
      break;
    case 2:
      /*lognormal (log10) case*/
      r = (double)sqrt(stat.variance[nblockv-1])/stat.mean[nblockm-1];
      r *= r;
      moy = (double)log10(stat.mean[nblockm-1]/sqrt(1.0+r));
      var = (double)log10(1.0+r);
      break;
    default:  
      Error("Type not defined in addstat");
      break;
    }

 

    (*realout).vector[i] = (double)sqrt(var)*(*realin).vector[i]+moy;
  }

  return;
}