Distribution de gauss

Signaler
Messages postés
3
Date d'inscription
samedi 15 novembre 2008
Statut
Membre
Dernière intervention
8 juin 2009
-
Messages postés
966
Date d'inscription
samedi 3 avril 2004
Statut
Membre
Dernière intervention
4 mars 2010
-
j'aimerai bien générer une distribution de rayons par la méthode de gauss et après avoir écris le code,  le resultat de la compilation affiche le fichier texte avec aucune valeur. Je ne sais pas pourquoi, svp quelqu'un pourrait-il m'aider?

 Voici le code joint:

#include <cstdlib>
#include <stdio.h>

#include <ctype.h>
#include <math.h>

#include <time.h>

#include "hmath.h"

/*  ****************Constant definition******************* */

const double mu = 4.93315;     /*mean value of the charge used
                           in the charge distribution function*/
const int nop = 5000;

const double sigma = 0.17352;
const double y00 = 5.271;
const double A = 20.83073;

/*************************************************************/
/* *****************  global variables ********************* */

double *charge;  /*charge distribution array*/
double *mass;     /*mass distribution array*/

double m0;  /*mean value of the mass*/
double Q0;  /*mean value of the charge*/

FILE *fp;

    
/************************** initialising particle mass and charge *************/
/* void init_part_mass_charge(double sigma)
{
     long int i;
    
    charge = new double [nop];
    mass = new double [nop];
       
          for (i=0;i<nop;i++){
              charge[i]= mu+sigma*sqrt(2.0)*
                        errorFn_inv(2.0*((double)rand()/(double)RAND_MAX)-1.0);
              mass[i] = charge[i]*charge[i]*charge[i];
             // printf("%ld %.3e %.3e\n", i, charge[i], mass[i]);
         }
}     */

void init_part_mass_charge(double sigma)
{
     long int i;
    
    charge = new double [nop];
    mass = new double [nop];
       
          for (i=0;i<nop;i++){
              charge[i]= mu+sigma/sqrt(2.0)*
                        errorFn_inv((2.0*y00*((double)rand()/(double)RAND_MAX)-y00));
              mass[i] = charge[i]*charge[i]*charge[i];
             // printf("%ld %.3e %.3e\n", i, charge[i], mass[i]);
         }
}    
/* *********************************************************************** */
/* ***************** MAIN ******************************************* */

int main(int argc, char* argv[]) {
      long int i;      
   
    srand(time(NULL));      /*permet d'initialiser le générateur de nombres
                            pseudo-aléatoires avec une graine différente*/
      
    init_part_mass_charge(sigma);
     
    /********* saving  data ******************/
   
    fp=fopen("distri.dat", "w");
   
    for (i = 0; i < nop ; i++){
    fprintf(fp, " %.3e %.3e\n",  charge[i], mass[i]);
    }
   
    fclose(fp);
    printf(" OK\n");
    delete [] charge;
    delete [] mass;
        
    system("PAUSE");
   
}

voici aussi le fichier math.h

  /* *************************************************** */
 /* ********************  hmath.h  ******************** */
/* **** My own mathematical function compendium ****** */

// Math. constant definitions:
const double   pi       = 3.14159265358979323846264338327950;

// Inverse error function
// Input  (double):    -1 <   y    < 1
// Output (double):  -inf < erf^-1 < inf
double errorFn_inv(double y) {
 
  double s, t, u, w, x, z;
  double k = y;     // store y before switching its sign
 
  if (y == 0)
    {
      x = 0;
      return x;
    }
  if (y > 1.0)
    {
      x = -log(0);     // to generate +inf
      return x;
    }
  if (y < -1.0)
    {
      x = log(0);     // to generate -inf
      return x;
    }
  if (y < 0)
    y = -y;     // switch the sign of y it's negative hence the comupation is for y >0
 
  z = 1.0 - y;
  w = 0.916461398268964 - log(z);
  u = sqrt(w);
  s = (log(u) + 0.488826640273108) / w;
  t = 1 / (u + 0.231729200323405);
  x = u * (1.0 - s * (s * 0.124610454613712 + 0.5)) -
    ((((-0.0728846765585675 * t + 0.269999308670029) * t +
       0.150689047360223) * t + 0.116065025341614) * t +
     0.499999303439796) * t;
  t = 3.97886080735226 / (x + 3.97886080735226);
  u = t - 0.5;
  s = (((((((((0.00112648096188977922 * u +
           1.05739299623423047e-4) * u - 0.00351287146129100025) * u -
         7.71708358954120939e-4) * u + 0.00685649426074558612) * u +
       0.00339721910367775861) * u - 0.011274916933250487) * u -
     0.0118598117047771104) * u + 0.0142961988697898018) * u +
       0.0346494207789099922) * u + 0.00220995927012179067;
  s = ((((((((((((s * u - 0.0743424357241784861) * u -
         0.105872177941595488) * u + 0.0147297938331485121) * u +
           0.316847638520135944) * u + 0.713657635868730364) * u +
         1.05375024970847138) * u + 1.21448730779995237) * u +
       1.16374581931560831) * u + 0.956464974744799006) * u +
     0.686265948274097816) * u + 0.434397492331430115) * u +
       0.244044510593190935) * t -
    z * exp(x * x - 0.120782237635245222);
  x += s * (x * s + 1.0);
 
  if(k < 0)
    return -x;
  else
    return x;
}
A voir également:

1 réponse

Messages postés
966
Date d'inscription
samedi 3 avril 2004
Statut
Membre
Dernière intervention
4 mars 2010
4
ben j'ai testé, et ça remplit bien le fichier (je précise, OS WinXP, compilateur VCE008)
essaie de débuger en vérifiant :
- que ton tableau est rempli correctement
- que ton fichier est ouvert correctement
...