
#include "calculQr.h"


/******************************************************************************/
/* VARIABLES GLOBALES                                                         */
/******************************************************************************/

double    dx;
long      ny;
double    *Qr;
forme_t   *formes;
ellipse_t *ellipses;
shaper_t  *shapers;


/******************************************************************************/
/* CALCUL DES ARETES D'UN POLYGONE                                            */
/*                                                                            */
/* Remplit la structure arete_t pour chaque arete d'un polygone.              */
/******************************************************************************/

void calcul_aretes(forme_t *forme,arete_t *aretes) {

  /* variables locales */
  
  long i,j;
  double *x,*y;


  /* vecteurs de points à utiliser */
  
  if (forme->x_rotation) {
    x=forme->x_rotation;
    y=forme->y_rotation;
  } else {
    x=forme->x;
    y=forme->y;
  }


  for (i=0;i<forme->nb;i++) {
    j=(i+1)%forme->nb;
    if (y[j]==y[i]) y[i]-=EPSILON;
    if (x[j]!=x[i]) {
      aretes[i].a=(y[j]-y[i])/(x[j]-x[i]);
    } else {
      aretes[i].a=INFINI;
    }
    aretes[i].b=y[i]-aretes[i].a*x[i];
    aretes[i].y_min=MIN(y[i],y[j]);
    aretes[i].y_max=MAX(y[i],y[j]);
  }
}


/******************************************************************************/
/* FONCTION DE MOTIF (CARRE, ANGLE FIXE)                                      */
/******************************************************************************/

void f_carre_fixe(position_t *position,tirage_t *tirage) {

  /* variables locales */

  double xi,yi,ri,wi;
  long i,j;


  xi=tirage->x;
  yi=tirage->y;
  ri=tirage->r;
  wi=tirage->w;
  for (i=position->dix;i<=position->fix;i++) {
    for (j=position->diy;j<=position->fiy;j++) {
      Qr[i*ny+j]*=wi;
    }
  }
}


/******************************************************************************/
/* FONCTION DE MOTIF (CARRE, ANGLE VARIABLE)                                  */
/******************************************************************************/

void f_carre_variable(position_t *position,tirage_t *tirage) {

  /* variables locales */

  double xi,yi,ri_2rc2,wi,cos_thetai,sin_thetai;
  long i,j;


  xi=tirage->x;
  yi=tirage->y;
  ri_2rc2=tirage->r*C_1_2RC2; /* ri/(2*sqrt(2)) */
  wi=tirage->w;
  cos_thetai=cos(tirage->theta); sin_thetai=sin(tirage->theta);
  for (i=position->dix;i<=position->fix;i++) {
    for (j=position->diy;j<=position->fiy;j++) {
      if (fabs(cos_thetai*(i*dx-xi)+sin_thetai*(j*dx-yi))<=ri_2rc2 &&
          fabs(sin_thetai*(i*dx-xi)-cos_thetai*(j*dx-yi))<=ri_2rc2) 
        Qr[i*ny+j]*=wi;
    }
  }
}


/******************************************************************************/
/* FONCTION DE MOTIF (CERCLE)                                                 */
/******************************************************************************/

void f_cercle(position_t *position,tirage_t *tirage) {

  /* variables locales */

  double xi,yi,ri2_4,wi;
  long i,j;


  xi=tirage->x;
  yi=tirage->y;
  ri2_4=.25*CARRE(tirage->r);
  wi=tirage->w;
  for (i=position->dix;i<=position->fix;i++) {
    for (j=position->diy;j<=position->fiy;j++) {
      if (CARRE(i*dx-xi)+CARRE(j*dx-yi)<ri2_4) Qr[i*ny+j]*=wi;
    }
  }
}


/******************************************************************************/
/* FONCTION DE MOTIF (SHAPER cos^2(pi.r) )                                    */
/******************************************************************************/

void f_shaper_cos_carre(position_t *position,tirage_t *tirage) {

  /* variables locales */

  double xi,yi,ri,ri2_4,wi,r2;
  long i,j;


  xi=tirage->x;
  yi=tirage->y;
  ri=tirage->r;
  ri2_4=.25*CARRE(ri);
  wi=tirage->w;
  for (i=position->dix;i<=position->fix;i++) {
    for (j=position->diy;j<=position->fiy;j++) {
      r2=CARRE(i*dx-xi)+CARRE(j*dx-yi);
      if (r2<ri2_4)
        Qr[i*ny+j]+=log((wi?wi:EPSILON))*CARRE(cos(T_PI*sqrt(r2)/ri));
    }
  }
}


/******************************************************************************/
/* FONCTION DE MOTIF (SHAPER exp[ 1 - 1/(1-(2r)^4) ])                         */
/******************************************************************************/

void f_shaper_exp(position_t *position,tirage_t *tirage) {

  /* variables locales */

  double xi,yi,ri,ri2_4,wi,r2;
  long i,j;


  xi=tirage->x;
  yi=tirage->y;
  ri=tirage->r;
  ri2_4=.25*CARRE(ri);
  wi=tirage->w;
  for (i=position->dix;i<=position->fix;i++) {
    for (j=position->diy;j<=position->fiy;j++) {
      r2=CARRE(i*dx-xi)+CARRE(j*dx-yi);
      if (r2<ri2_4)
        Qr[i*ny+j]+=log((wi?wi:EPSILON))*exp(1.-1./(1.+EPSILON-CARRE(r2/ri2_4)));
    }
  }
}


/******************************************************************************/
/* FONCTION DE MOTIF (POLYGONE)                                               */
/*                                                                            */
/* Un point de la matrice Qr est dans un polygone si le nombre d'intersections*/
/* entre une demi-droite partant de ce point et les aretes du polygone est im-*/
/* pair.                                                                      */
/******************************************************************************/

void f_polygone(position_t *position,tirage_t *tirage) {

  /* variables locales */

  double xi,yi,ri,wi,cos_thetai,sin_thetai,x,y,yc;
  long i,j,k,nb;
  forme_t forme=formes[tirage->numero_motif];
  arete_t *aretes=forme.intervalles.aretes;


  xi=tirage->x;
  yi=tirage->y;
  ri=tirage->r;
  wi=tirage->w;


  /* calculs des segments du polygone si rotation */

  if (forme.x_rotation) {
    cos_thetai=cos(tirage->theta); sin_thetai=sin(tirage->theta);
    for (i=0;i<forme.nb;i++) {
      forme.x_rotation[i]= forme.x[i]*cos_thetai+forme.y[i]*sin_thetai;
      forme.y_rotation[i]=-forme.x[i]*sin_thetai+forme.y[i]*cos_thetai;
    }
    calcul_aretes(&forme,aretes);
  }


  for (i=position->dix;i<=position->fix;i++) {
    for (j=position->diy;j<=position->fiy;j++) {
    
      /* nombre de segments intersectes */

      nb=0;
      x=(i*dx-xi);
      y=(j*dx-yi);
      for (k=0;k<forme.nb;k++) {
        yc=aretes[k].a*x+aretes[k].b*ri;
        if (y<=yc && yc>=ri*aretes[k].y_min
                  && yc<=ri*aretes[k].y_max) {
          nb++;
        }
      }


      /* si le point est dans le polygone -> mise à jour de Qr */

      if ((nb%2)==1) Qr[i*ny+j]*=wi;
    }
  }
}


/******************************************************************************/
/* FONCTION DE MOTIF (SPLINE)                                                 */
/*                                                                            */
/* Un point de la matrice Qr est dans une spline fermee si le nombre d'inter- */
/* sections entre une demi-droite partant de ce point et la spline est impair.*/
/******************************************************************************/

void f_spline(position_t *position,tirage_t *tirage) {

  /* variables locales */

  double xi,yi,ri,wi,cos_thetai,sin_thetai,*sx,*sy,x,y,yc,r[3];
  long i,j,k,l,nb,nb_racines;
  forme_t forme=formes[tirage->numero_motif];
  polynome_t py,*polynomes_x=forme.intervalles.spline.x,
                *polynomes_y=forme.intervalles.spline.y;


  xi=tirage->x;
  yi=tirage->y;
  ri=tirage->r;
  wi=tirage->w;


  /* calculs des coefficients des polynomes si rotation */

  if (forme.x_rotation) {
    cos_thetai=cos(tirage->theta); sin_thetai=sin(tirage->theta);
    for (i=0;i<forme.nb;i++) {
      forme.x_rotation[i]= forme.x[i]*cos_thetai+forme.y[i]*sin_thetai;
      forme.y_rotation[i]=-forme.x[i]*sin_thetai+forme.y[i]*cos_thetai;
    }
    coeffs_spline(forme.x_rotation,polynomes_x,forme.nb);
    coeffs_spline(forme.y_rotation,polynomes_y,forme.nb);
    sx=forme.x_rotation;
    sy=forme.y_rotation;
  } else {
    sx=forme.x;
    sy=forme.y;
  }


  for (i=position->dix;i<=position->fix;i++) {
    for (j=position->diy;j<=position->fiy;j++) {

      /* nombre d'intersections avec la courbe */

      nb=0;
      x=(i*dx-xi);
      y=(j*dx-yi);
      for (k=0;k<forme.nb-1;k++) {
        nb_racines=racines(sx[k]-x/ri,polynomes_x[k],r);
        py=polynomes_y[k];
        for (l=0;l<nb_racines;l++) {
          if (r[l]>=0. && r[l]<=1.) {
            yc=ri*(sy[k]+r[l]*(py.c1+r[l]*(py.c2+r[l]*py.c3)));
            if (y<=yc) nb++;
          }
        }
      }


      /* si le point est dans la spline -> mise à jour de Qr */

      if ((nb%2)==1) Qr[i*ny+j]*=wi;

    }
  }
}


/******************************************************************************/
/* FONCTION DE MOTIF (ELLIPSE)                                                */
/******************************************************************************/

void f_ellipse(position_t *position,tirage_t *tirage) {

  /* variables locales */

  double xi,yi,ri,wi,angle,c_ri_cos_phi,c_ri_sin_phi,ri_grand_axe,x,y;
  long i,j;
  ellipse_t ellipse=ellipses[tirage->numero_motif];


  xi=tirage->x;
  yi=tirage->y;
  ri=tirage->r;
  wi=tirage->w;
  angle=(tirage->angle==fixe?ellipse.angle:tirage->theta);
  c_ri_cos_phi=ellipse.foyer*ri*cos(angle);
  c_ri_sin_phi=ellipse.foyer*ri*sin(angle);
  ri_grand_axe=2.*ri*ellipse.demi_grand_axe;
  for (i=position->dix;i<=position->fix;i++) {
    for (j=position->diy;j<=position->fiy;j++) {
      x=i*dx-xi;
      y=j*dx-yi;
      if (sqrt(CARRE(x-c_ri_cos_phi)+CARRE(y-c_ri_sin_phi))
         +sqrt(CARRE(x+c_ri_cos_phi)+CARRE(y+c_ri_sin_phi))<=ri_grand_axe)
        Qr[i*ny+j]*=wi;
    }
  }
}


/******************************************************************************/
/* FONCTION DE MOTIF (SHAPER QUELCONQUE)                                      */
/******************************************************************************/

void f_shaper_quelconque(position_t *position,tirage_t *tirage) {

  /* variables locales */

  double xi,yi,ri,ri2_4,wi,r2,r,x;
  long i,j,inf;
  shaper_t shaper=shapers[tirage->numero_motif];


  xi=tirage->x;
  yi=tirage->y;
  ri=tirage->r;
  ri2_4=.25*CARRE(ri);
  wi=tirage->w;
  for (i=position->dix;i<=position->fix;i++) {
    for (j=position->diy;j<=position->fiy;j++) {
      r2=CARRE(i*dx-xi)+CARRE(j*dx-yi);
      if (r2<ri2_4) {
        r=sqrt(r2)/ri;
        x=2.*r*(shaper.nb-1);
        inf=(long)x;
        Qr[i*ny+j]+=log((wi?wi:EPSILON))
                   *(shaper.points[inf]+(x-inf)*(shaper.points[inf+1]-shaper.points[inf]));
      }
    }
  }
}


/******************************************************************************/
/* CALCUL DE QR                                                               */
/*                                                                            */
/* Construction de la matrice Qr en effectuant pour chacun de ses points le   */
/* produit des multiplieurs des points generes aleatoirement dont le motif    */
/* 'eclaire' le point. On commence par parcourir les motifs de type fonction  */
/* lissante.                                                                  */
/******************************************************************************/

void calcul_Qr(double _dx,long nx,long _ny,double *_Qr,forme_t *_formes,
               ellipse_t *_ellipses,shaper_t *_shapers,long debut_non_shaper,
               positions_t positions,tirage_t *tirages) {

  /* variables locales */

  long i;
  void (*f_motif[9])(position_t *,tirage_t *)={
    f_carre_fixe,f_carre_variable,f_cercle,f_shaper_cos_carre,f_shaper_exp,
    f_polygone,f_spline,f_ellipse,f_shaper_quelconque
  };
  

  /* initialisation des variables globales */

  dx=_dx;
  ny=_ny;
  Qr=_Qr;
  formes=_formes;
  ellipses=_ellipses;
  shapers=_shapers;


  /* fonctions de motif de type fonction lissante */

  if (debut_non_shaper>0) {
    for (i=0;i<nx*ny;i++) Qr[i]=0.;
    for (i=0;i<debut_non_shaper;i++)
      f_motif[tirages[i].numero_fonction](positions.tableau+i,tirages+i);
    for (i=0;i<nx*ny;i++) Qr[i]=exp(Qr[i]);
  } else for (i=0;i<nx*ny;i++) Qr[i]=1.;


  /* autres fonctions de motif */

  for (i=debut_non_shaper;i<positions.n;i++)
    f_motif[tirages[i].numero_fonction](positions.tableau+i,tirages+i);
}

