
#include "tirages.h"


/******************************************************************************/
/* VALEUR TIREE SELON UNE LOI NORMALE CENTREE REDUITE                         */
/******************************************************************************/

double randn(void) {

  /* variables locales */

  double x1,x2,w;
  static unsigned char valeur_disponible=0;
  static double autre_valeur;


  if (!valeur_disponible) {

    /* tirage de deux valeurs suivant une loi normale */
    /* centree reduite et retour d'une d'elle         */

    do {
      x1=2.*randu()-1.;
      x2=2.*randu()-1.;
      w=CARRE(x1)+CARRE(x2);
    } while (w>=1.);
    w=sqrt((-2.*log(w))/w);
    autre_valeur=x1*w;
    valeur_disponible=1;
    return x2*w;

  } else {

    /* retour d'une valeur precedemment calculee */

    valeur_disponible=0;
    return autre_valeur;

  }
}


/******************************************************************************/
/* ZONES D'INFLUENCE                                                          */
/*                                                                            */
/* Determination des points de la matrice Qr contenus dans un carre de centre */
/* (xi,yi) et de cote ri (l'angle du carre par rapport a Qr etant nul).       */
/******************************************************************************/

long zone_influence(double xi,double yi,double ri,long nx,long ny,
                    double dx,positions_t *positions) {

  /* variables locales */

  long dix,fix,diy,fiy,retour=0;


  /* limites de la zone d'influence */

  dix=MAX(   0,(long)ceil ((xi-ri/2.)/dx));
  fix=MIN(nx-1,(long)floor((xi+ri/2.)/dx));


  /* si le point a une zone d'inflence selon x, on continue */

  if (dix<=fix) {

    diy=MAX(   0,(long)ceil ((yi-ri/2.)/dx));
    fiy=MIN(ny-1,(long)floor((yi+ri/2.)/dx));


    /* si le point a une zone d'inflence selon y, on continue */

    if (diy<=fiy) {

      /* sauvegarde des donnees */

      positions->tableau[positions->n].x=xi;
      positions->tableau[positions->n].y=yi;
      positions->tableau[positions->n].r=ri;
      positions->tableau[positions->n].dix=dix;
      positions->tableau[positions->n].fix=fix;
      positions->tableau[positions->n].diy=diy;
      positions->tableau[positions->n].fiy=fiy;
      positions->n++;
      retour=1;

    }
  }
  return retour;
}


/******************************************************************************/
/* TIRAGE DES POSITIONS                                                       */
/*                                                                            */
/* Generation d'un point (xi,yi,ri) et determination de sa zone d'influence.  */
/******************************************************************************/

void tirer_positions(long N,double x_min,double x_max,double y_min,double y_max,
                     double u_max,double r_min,double bord,double dx,
                     positions_t *positions,long nx,long ny) {

  /* variables locales */

  long i;
  double xi,yi,ui,ri,delta_x=x_max-x_min,delta_y=y_max-y_min;


  /* tirage des positions */

  for (i=0;i<N;i++) {

    /* position */

    xi=x_min-.5*bord+(delta_x+bord)*randu();
    yi=y_min-.5*bord+(delta_y+bord)*randu();
    ui=                       u_max*randu();
    ri=1./sqrt(1./CARRE(r_min)-ui);


    /* calcul de la zone d'influence et sauvegarde si elle est valide */

    zone_influence(xi,yi,ri,nx,ny,dx,positions);

  }
}


/******************************************************************************/
/* TIRAGE DES MULTIPLIEURS                                                    */
/*                                                                            */
/* Generation d'un multiplieur pour chaque point (xi,yi,ri) en fonction de la */
/* loi, ainsi que d'un angle de rotation pour le motif. On construit la liste */
/* des variables aleatoires retournees sous Matlab.                           */
/******************************************************************************/

void tirer_multiplieurs(positions_t positions,tirage_t *tirages,loi_t loi,
                        motif_t type_motif,angle_t angle,numero_fonction_t numero_fonction,
                        long numero_motif,long numero_zone) {

  /* variables locales */

  long i;
  double cste;
  parametres_t p=loi.parametres;


  /* recopie des nouvelles positions, tirage des angles de   */
  /* rotation et affectation des caracteristiques de la zone */

  for (i=0;i<positions.n;i++) {
    memcpy(&(tirages[i]),&(positions.tableau[i]),3*sizeof(double));
    tirages[i].theta=2.*T_PI*randu();
    tirages[i].type_motif     =type_motif;
    tirages[i].angle          =angle;
    tirages[i].numero_fonction=numero_fonction;
    tirages[i].numero_motif   =numero_motif;
    tirages[i].numero_zone    =numero_zone;
  }


  /* tirage des multiplieurs */
  
  switch (loi.type) {
    case log_normal:
      cste=sqrt(p.ln.sigma2);
      for (i=0;i<positions.n;i++) tirages[i].w=exp(cste*randn()+p.ln.mu);
      break;
    case log_exponential:
      cste=pow(1.+p.le.Temp,1./p.le.Temp);
      for (i=0;i<positions.n;i++) tirages[i].w=pow(cste*randu(),p.le.Temp);
      break;
    case log_poisson:
      for (i=0;i<positions.n;i++) tirages[i].w=p.lp.W0;
      break;
  }
}

