source: Sophya/trunk/FrEROS/AnaLC/nbgene.cc@ 3891

#include "machdefs.h"
#include "sopnamsp.h"  // using SOPHYA namespace
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "fsvcache.h"
#include "nbtri.h"
#include "nbmath.h"
#include "nbinteg.h"
#include "strutil.h"
#include "fsvst.h"
#include "srandgen.h"
#include "nbsread.h"
#include "nbgene.h"

/*==========================================================================*/
int Give_Annee_Obs(int lp)
/* retourne l'annee d'observation codee en entier du fichier de suivi utilise */
{
int j,anneeobs;
double t=0.;
anneeobs = -1;
for(j=0;j<nmes[0];j++) {
  t = timeu[0][j].TStart / JourSec;
  if( t <= 0. ) continue;
  /* confusion possible smc 93 avec debut lmc 9394 */
  if( t > 1288. && t < 1314. ) continue;
  if( t >  700. && t < 850.  ) {anneeobs=9192; break;}
  if( t >  950. && t < 1200. ) {anneeobs=9293; break;}
  if( t > 1300. && t < 1570. ) {anneeobs=9394; break;}
  if( t > 1630. && t < 1900. ) {anneeobs=9495; break;}
  if( t > 2300. && t < 9999. ) {anneeobs=9697; break;}
}
if(lp>0) printf("annee_obs = %d  (j=%d t=%f)\n",anneeobs,j,t);
return (anneeobs);
}

/*=========================================================================*/
float Rayon_SP_Th(float magv)
/* Fournit le rayon de l'etoile de SP en rayons solaires
*  en fonction de sa magnitude Vj.
*  Formule theorique uniquement pour la sequence principale:
*  Rstar_SP = exp(-0.22 * (magv - MagV_sol) ) avec MagV_sol = 5.1
*/
{
return( (float) exp(-0.22*(magv-5.1)) );
} 

/*==========================================================================*/
void AMPLML(int io)
/* Generation MonteCarlo */
{
int i,ic,j,tf;
float x;
double t,tt,u2,u;
int Nseq_tf = 3;
unsigned short seed_16v[3];

mc.TauSim = mc.TauSimI;
mc.T0Sim = mc.T0SimI;
mc.U0Sim = mc.U0SimI;
mc.USim = mc.USimI;

if(io == 1) {

  /* initialisation des generateurs */
  seed_16v[0]=(unsigned short) mc.iseed1;
  seed_16v[1]=(unsigned short) mc.iseed2;
  seed_16v[2]=(unsigned short) mc.iseed3;
  RandGen->SetSeed(seed_16v);

  mc.NGene = mc.NPoints_ampli = mc.NPoints_ampli_tf = mc.Net_ampli_tf = 0;

  if(mc.montecar<=0)   return;
  if ( mc.Taille_Finie > 2 )
    { printf("mauvaise valeur mc.Taille_Finie %d\n",mc.Taille_Finie);
    exit(-1); }
  if ( mc.Blending > 2 || mc.Blending < 0 )
    { printf("mauvaise valeur mc.Blending %d\n",mc.Blending);
    exit(-1); }

  printf("MONTECARLO %2d: seed %d %d %d\n",io,mc.iseed1,mc.iseed2,mc.iseed3);
  printf("   Generation %d, Taille_Finie %d blending %d\n"
        ,mc.montecar,mc.Taille_Finie,mc.Blending);
  printf("   Nombre de generations par etoile %d\n",mc.NGenStar);
  printf("   TauSim=%f TauMinI=%f TauMaxI=%f TireTau=%p\n"
        ,mc.TauSimI,mc.TauMinI,mc.TauMaxI,mc.TireTau);
  printf("   T0Sim=%f T0MinI=%f T0MaxI=%f TireT0=%p\n"
        ,mc.T0SimI,mc.T0MinI,mc.T0MaxI,mc.TireT0);
  printf("   U0Sim=%f U0MinI=%f U0MaxI=%f TireU0=%p\n"
        ,mc.U0SimI,mc.U0MinI,mc.U0MaxI,mc.TireU0);
  printf("   USim =%f UMinI =%f UMaxI =%f TireU=%p\n"
        ,mc.USimI,mc.UMinI,mc.UMaxI,mc.TireU);
  printf("   coupure pour taille finie a U/u0 < %f\n",mc.Usu0_Cut);

  if( (mc.U0MaxI <= mc.U0MinI && mc.U0SimI <= 0. &&  mc.Taille_Finie==0 )
   || (mc.U0MaxI <= mc.U0MinI && mc.U0SimI < 0.  &&  mc.Taille_Finie>0  ) ) {
    printf("Mauvaises valeurs U0SimI U0MinI U0MaxI\n");
    exit(-1);
  }
  if( mc.UMaxI <= mc.UMinI && mc.USimI <= 0. && mc.Taille_Finie>0 ) {
    printf("Mauvaises valeurs USimI UMinI UMaxI\n");
    exit(-1);
  }
  if( mc.TauMaxI <= mc.TauMinI && mc.TauSimI <= 0.) {
    printf("Mauvaises valeurs TauSimI TauMinI TauMaxI\n");
    exit(-1);
  }

} else if(io == 2) {

  if( mc.T0MaxI <= mc.T0MinI && mc.T0SimI <= 0. ) {
    mc.T0MinI_eff=TFIRST;
    mc.T0MaxI_eff=TLAST;
  } else {
    mc.T0MinI_eff=mc.T0MinI;
    mc.T0MaxI_eff=mc.T0MaxI;
  }

  if(mc.montecar<=0) return;

  if( mc.Blending>0 && mc.Get_Par_Blending==0 ) {
    printf("Vous avez demande une generation avec du Blending\n");
    printf(" mais vous n'avez pas fourni la generation des parametres\n");
    exit(-1);
  }

  printf("MONTECARLO %2d: ccd %6d, T0MinI_eff=%f T0MaxI_eff=%f\n"
         ,io,ccdnum,mc.T0MinI_eff,mc.T0MaxI_eff);

} else if(io == 3) {

  mc.NGene++;

  /* Blending preparation */
  if( mc.Blending>0 ) mc.Get_Par_Blending();

  /* Tirage aleatoire de u0 */
  if ( mc.TireU0 != 0 ) {
      x = drand01();
      mc.U0Sim = mc.TireU0(x);
  } else if(mc.U0MaxI > mc.U0MinI) {
    mc.U0Sim= -1.;
    while( mc.U0Sim <= 0. ) {
      mc.U0Sim = mc.U0MinI + (mc.U0MaxI-mc.U0MinI)*drand01();
      if(mc.U0Sim<=0.) printf("ATTENTION: U0 SIMULE A %f min=%f max=%f\n"
                             ,mc.U0Sim,mc.U0MinI,mc.U0MaxI);
  } }

  /* generation plate de tau */
  if ( mc.TireTau != 0 ) {
     x = drand01();
     mc.TauSim = mc.TireTau(x);
  } else if ( mc.TauMaxI > mc.TauMinI ) {
     mc.TauSim = mc.TauMinI
               + (mc.TauMaxI-mc.TauMinI)*drand01();
  }

  /* Tirage aleatoire de t0 */
  if ( mc.TireT0 != 0 ) {
     x = drand01();
     mc.T0Sim = mc.TireT0(x);
  } else if( mc.T0MaxI_eff > mc.T0MinI_eff ) {
     mc.T0Sim = mc.T0MinI_eff
              + (mc.T0MaxI_eff-mc.T0MinI_eff)*drand01();
  }

  /* generation plate de U */
  if( mc.Taille_Finie > 0 ) {
    if ( mc.TireU != 0 ) {
       x = drand01();
       mc.USim = mc.TireU(x);
    } else if ( mc.Taille_Finie>0 && mc.UMaxI > mc.UMinI ) {
       mc.USim= mc.UMinI + (mc.UMaxI-mc.UMinI)*drand01();
    }
    if(mc.Taille_Finie==2) {
      if( mc.Blending==0 ) mc.USim *= starcal.Rstar;
      else if( mc.Blending==1 ) mc.USim *= mc.starcal1.Rstar;
      else if( mc.Blending==2 ) mc.USim *= mc.starcal2.Rstar;
    }
  }

  /* calcul amplification ponctuelle/taille finie (paczinsky) */
  mc.A0MaxPt = GRAND;
  if( mc.U0Sim > 0. ) mc.A0MaxPt = ampli_ponct(mc.U0Sim);
  mc.A0Max = (mc.Taille_Finie && mc.USim>0.) ?
             ampli_tf(mc.U0Sim,mc.USim,Nseq_tf) : mc.A0MaxPt;
  if( mc.Blending>0) for(ic=0;ic<NCOULMX;ic++) {
    if( mc.Blending == 1 )
        mc.A0Max_blend[ic] = (mc.A0Max -1.)* mc.coeff_Arec1[ic] + 1.;
    else if( mc.Blending == 2 )
        mc.A0Max_blend[ic] = (mc.A0Max -1.)* mc.coeff_Arec2[ic] + 1.;
  }

  tf = 0;
  for(ic=0;ic<NCOULMX;ic++) {
    for(i=0;i<nmes[ic];i++) {
      t = date[ic][i];
      tt = (t-mc.T0Sim)/mc.TauSim;
      u2 = mc.U0Sim*mc.U0Sim + tt*tt;
      u = sqrt(u2);
      ampli[ic][i] = ampli_ponct(u);
      mc.NPoints_ampli++;
      if( mc.Taille_Finie && t<GRAND2 && mc.USim>0. ) {
        if( mc.USim > mc.Usu0_Cut * u ) {
          ampli[ic][i] = ampli_tf(u,mc.USim,Nseq_tf);
          mc.NPoints_ampli_tf++;
          if(tf==0) { mc.Net_ampli_tf++ ; tf = 1;}
        }
      }
      if(mc.Blending==1)
             ampli[ic][i] = (ampli[ic][i] -1.)* mc.coeff_Arec1[ic] + 1.;
        else if(mc.Blending==2)
             ampli[ic][i] = (ampli[ic][i] -1.)* mc.coeff_Arec2[ic] + 1.;
  } }

  if(debug>1) {
    for(ic=0;ic<NCOULMX;ic++) {
      printf("ampli[%3d]-1 u0=%12.5g t0=%12.5g tau=%12.5g U=%12.5g:\n"
            ,ic,mc.U0Sim,mc.T0Sim,mc.TauSim,mc.USim);
      for(i=0;i<nmesure[ic];i++)
        { j=indexu[ic][i]; printf(" %12.5g",ampli[ic][j]-1.);
        if(i%10==9) printf("\n");}
      if((nmesure[ic]-1)%10!=9) printf("\n");
  } }

} else if(io == 4) {

  return;

} else if(io == 5) {

  if(mc.montecar<=0) return;

  printf("MONTECARLO %2d: Nombre d'evenements engendres %d\n",io,mc.NGene);
  printf("Nb de pts amplifies %d, nb de pts amplifies avec taille finie %d\n"
        ,mc.NPoints_ampli,mc.NPoints_ampli_tf);
  printf("Nb d'etoiles amplifiees avec taille finie %d\n"
        ,mc.Net_ampli_tf);
  RandGen->SetSeed(seed_16v);
  mc.iseed1 = (int_4)seed_16v[0];
  mc.iseed2 = (int_4)seed_16v[1];
  mc.iseed3 = (int_4)seed_16v[2];
  printf("Aleatoire_de_fin_de_Generation %d %d %d\n"
        ,mc.iseed1,mc.iseed2,mc.iseed3);

}
return;
}

/*==========================================================================*/
double ampli_ponct(double u)
/* amplification paczinsky dans le cas d'une source ponctuelle */
{
u *= u;
if(u==0.) return(GRAND2);
u = (u+2.)/sqrt(u*(u+4.));
return (u);
}

/*==========================================================================*/
double ampli_tf(double u, double U,int Nseq)
/*
Amplification paczinsky dans le cas d'une source non ponctuelle.
  ( integration analytique partie interieure et
    integration numerique partie exterieure )
Integration exterieure par methode de Gauss (Higher order terms).
  u = parametre d'impact en rayon d'Einstein en fct du temps
  U = rayon de l'etoile projete en rayon d'Einstein
  Nseq = nombre de pas pour l'integrale
*/
{
int i,j;
double u2,U2,lim1,lim2,dlim,A0,xc,x,nx,b,teta,Ka,z;
int mIhoqN;
double *IhoqNX, *IhoqNW;

if( Nseq <= 0 ) Nseq = 2;
if( u<0. || U<0. ) return(-1.);
if( u==0. && U==0. ) return(GRAND2);

mIhoqN=mIhoq9; IhoqNX=Ihoq9X; IhoqNW=Ihoq9W;

u2 = u*u;
U2 = U*U;

Ka = 0.;
lim1 = -1.;
lim2 = 1.;
dlim = (lim2 - lim1)/Nseq;

/*** CAS 1 : u<=U (u=0 possible si U!=0) ***                */
/* A = A0 + A1                                              */
/* A0 = partie integrable (sum(dx,0<x<U-u))                 */
/*    = (U-u) * sqrt((U-u)**2 + 4) / U2                     */
/* A1 = 2*u/(Pi*U2) * sum(dx,-1<x<1)                        */
/*    {(u2+2)/sqrt(u2+4)*acos[(u*(1+x**2)/2+U*x)/(U+x*u)]}  */

if(u<=U) {
  /* on calcule d'abord la partie analytique A0 */
  x = (U-u);
  A0 = x * sqrt(x*x + 4.);

  /* on integre ensuite la partie numerique A1 */
  for(j=0;j<Nseq;j++) {
    xc = lim1 + (j+0.5)*dlim;
    for(i=0;i<mIhoqN;i++) {
      x = xc + IhoqNX[i]*dlim;
      nx = U + x*u;
      teta = acos( ( u*(1.+x*x)/2. + U*x ) /nx );
      nx *= nx;
      b = (nx + 2.) / sqrt(nx + 4.);
      Ka += b * teta * IhoqNW[i];
  } }
  A0 += 2. * u / Pi * Ka * dlim;
  A0 /= U2;
}

/*** CAS 2 : u>U   -->  U peut tendre vers 0 ***               */
/* A = 2/(Pi*U) * sum(dx,-1<x<1)                               */
/*    {(u2+2)/sqrt(u2+4)*acos[1-U2*(1-x**2)/(2*u*(u+x*U))]}    */
/* quand U->0, la divergence (1/U) est enlevee en developpant  */
/* acos(1-z) = sqrt(2*z)*[1+z/12+3*z2/160+...]                 */

else if(u>U) {
  for(j=0;j<Nseq;j++) {
    xc = lim1 + (j+0.5)*dlim;
    for(i=0;i<mIhoqN;i++) {
      x = xc + IhoqNX[i]*dlim;
      nx = u + x*U;
      z = (1.-x*x) /(2.*u*nx);
      nx *= nx;
      b = (nx + 2.) / sqrt(nx + 4.);
      if(U > 0.001*u) {
        teta = acos( 1. - U2*z ) / U;
      } else {        /* cas U-> 0 */
        teta = sqrt(2.*z);
        z *= U2;
        teta *= 1. + z*(0.08333333333333 + 0.01875*z);
      }
      Ka += b * teta * IhoqNW[i];
  } }
  A0 = 2. / Pi * Ka * dlim;
}

/*** CAS IMPOSSIBLE ****/
else {
  printf("ampli_tf: CAS IMPOSSIBLE u=%f U=%f\n",u,U);
  A0 = -1.;
}

return (A0);
}

/*==========================================================================*/
double u0_ponct(double A)
/*
Parametre d'impact pour une amplification A
  dans le cas d'une source ponctuelle.
input: A = amplification
*/
{
double u0;
if(A<1.) return(-1.);
if(A==1.) return(GRAND2);
u0 = 2. * ( A/sqrt(A*A-1.) - 1.);
u0 = sqrt(u0);
return(u0);
}

/*==========================================================================*/
double u0_tf(double A,double U,double prec,int *niter,double *Afinal)
/*
Parametre d'impact pour une amplification A
  dans le cas d'une source non ponctuelle.
input:
  A = amplification
  U = Rstar projete / Reinst
  prec = precision souhaitee sur u0
  niter = nombre maximum d'iterations permises
output:
  niter = nombre d'iterations utilisees (-1 si pas converge)
  Afinal = amplification finale calculee pour u0
  return = u0 si trouve
           -1. si A<1.
           -2. si A=1.
           -3. si U<=0
           -4. si A trop grand
           -5. si encadrement de A non trouve
*/
{
int nitermx;
double a0=0,a1,a2,u0=0,u1,u2,test;

nitermx = *niter;
if(nitermx<=0) nitermx = 25;

*Afinal = A;
*niter = 0;
if(A<1.)  return(-1.);
if(A==1.) return(-2.);
if(U<=0.) return(-3.);

u1 = 0.;
a1 = ampli_tf(u1,U,-1);
if(a1==A) return(0.);
/* on ne peut atteindre la valeur A */
if(a1<A) { *Afinal = a1; return(-4.);}

/* on recherche un intervalle encadrant A */
u2 = u0_ponct(A);
a2 = ampli_tf(u2,U,-1);
while( a2>A && *niter<=nitermx ) {
  /* printf("recherche inter %f (%f), %f (%f)\n",u1,a1,u2,a2); */
  u1 = u2;
  a1 = a2;
  u2 *= 5.;
  a2 = ampli_tf(u2,U,-1);
  (*niter)++;
}

if( *niter>nitermx ) {
  *Afinal = a2;
  return(-5.);
}
  
/* a ce niveau u1 (a1) et u2 (a2) encadrent u0 (a0) */
test = 100.*prec;
*niter = 0;
while (test > prec && *niter<= nitermx ) {
  /* printf("recherche conv %f (%f), %f (%f)\n",u1,a1,u2,a2); */
  u0 = (u1+u2)/2.;
  a0 = ampli_tf(u0,U,-1);
  (*niter)++;
  if(a0<A) {
    u2 = u0;
    a2 = a0;
  } else {
    u1 = u0;
    a1 = a0;
  }
  test = fabs(a0-A);
}

if( *niter > nitermx ) *niter = -1;
/* printf("conv atteinte %f  (%f), %f (%f)\n",u1,u2,a2); */
*Afinal = a0;
return(u0);

}