//  ***********************************************
//  *  Corso Statistica L.Ramello - Esercizio 1   *
//  *                                             *
//  *  Generazione di eventi secondo una legge    *
//  *  di decadimento radioattivo.                *
//  *                                             *
//  *  Parametri:                                 *
//  *   alpha   - costante di decadimento         *
//  *   delta_t - intervallo di tempo di misura   *
//  *   N0      - numero iniziale di nuclei       *
//  *                                             *
//  *  Output:                                    *
//  *   N(t)    - nuclei sopravvissuti vs. t      *
//  *             (ntupla 1 in decrad1.hbook)     *
//  *                                             *
//  *   Per visualizzare il risultato:            *
//  *   PAW > nt/pl 1.Nuclei%Tempo                *
//  *         -> numero di nuclei sopravvissuti)  *
//  *   PAW > nt/pl 1.Num%Tempo                   *
//  *         -> ultimo numero casuale estratto   *
//  *            per ogni intervallo di tempo     *
//  *                                             *
//  *  M.Obertino, S.Valchierotti, A.Stamerra     *
//  *  8/9/1999                                   *
//  ***********************************************

#include "cfortran.h"
#include "packlib.h"
#include "kernlib.h"
#include "stdlib.h"
#define debug

#define PAWC_SIZE 50000

typedef struct { float PAW[PAWC_SIZE]; } PAWC_DEF;
#define PAWC COMMON_BLOCK(PAWC,pawc)
COMMON_BLOCK_DEF(PAWC_DEF,PAWC);
PAWC_DEF PAWC;

main()
{
  int hid=1,istat=0,icycle=0;
  int nvar=3;
  char Title[3][8]={"Tempo","Nuclei","Num"};
  float VecNtu[3];
  int record_size=1024;

  float vec[100];
  int len=100;
  int N0=100,count,t,i,N,j;
  float alpha=0.01,delta_t=1.;

  // Ntuple stuff

  HLIMIT(PAWC_SIZE);
  HROPEN(1,"decrad1","decrad1.hbook","N",record_size,istat);
  HBOOKN(hid,"Decadimento Radioattivo",nvar," ",5000,Title);

  printf("inizio danze ! \n");
  N=N0;
  for (t=0;t<300;t++) {  //--Loop sul tempo totale di osservazione (300 sec)

    count=0;

// Genera sequenza numeri casuali con funzione built-in
//    for (j=0;j<N;j++) {
//      vec[j]=rand(1.)/32768.;
//      VecNtu[2]=vec[j];
//    }

//  la chiamata RANMAR(vec,len) da' un Warning di compilazione ...
    RANMAR(vec[0],len);  //--Genera sequenza numeri casuali
    VecNtu[2]=vec[N-1];  //--Salva ultimo numero casuale utilizzato 

    for (i=0;i<N;i++) {  //--Loop sui nuclei sopravvissuti

      if(vec[i]<=(alpha*delta_t)) count=count+1;
    }
      N=N-count;
      VecNtu[0]=t;
      VecNtu[1]=N;

      HFN(hid,VecNtu);  //--Riempie t-esima riga della Ntupla
  }

  printf("fine danze! riempio file--> \ndecrad1.hbook\n");

  HROUT(0,icycle," ");
  HREND("decrad1");
  KUCLOS(1," ",1);

}