Sujet Précédent Sujet Suivant

Distribution de gauss

Flopy21 Messages postés 3 Date d'inscription samedi 15 novembre 2008 Statut Membre Dernière intervention 8 juin 2009 - 8 juin 2009 à 17:21
cs_juju12 Messages postés 966 Date d'inscription samedi 3 avril 2004 Statut Membre Dernière intervention 4 mars 2010 - 8 juin 2009 à 18:45
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;
}

1 réponse

Réponse 1 / 1
cs_juju12 Messages postés 966 Date d'inscription samedi 3 avril 2004 Statut Membre Dernière intervention 4 mars 2010 4
8 juin 2009 à 18:45
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
...
0

Discussions similaires

transposer une formule excel en vba?
cs_Fabian123 -
us_30 -

7 réponses
matrice et gauss
cs_tibo76 -
751 -

2 réponses
pivot de Gauss en Pascal!
bidibul12 -
bidibul12 -

2 réponses
algorithme de gauss et decomposition LU
speedamine -
ipower -

4 réponses
Echelle graphique normalisée
Herve_be -
BunoCS -

2 réponses
Rejoignez-nous