phasemaker.c 3.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #define pi 3.14159265358979323844
  5. #define NB 1
  6. void gauss(double *g1,double *g2);
  7. long IDUM=-1234;double ran2(void);
  8. #define NTYPES 5
  9. int main(){
  10. double p[4],r[4],mass[NTYPES]={938.3,939.6,139.58,139.58,493.7};
  11. int ident[NTYPES]={2212,2112,211,-211,321};
  12. double dummy,weight,pmag2,pmag,temp,rx,ry,rz,tau,cthet,sthet,phi;
  13. int i,id,alpha,iblankline,ib;
  14. char filename[80];
  15. FILE *fptr;
  16. printf("idents are %d %d %d\n",ident[0],ident[1],ident[2]);
  17. printf("What is the temperature? (in MeV)\n");
  18. scanf("%lf",&temp);
  19. printf("What are rx,ry,rz and tau? (in fm)\n");
  20. scanf("%lf %lf %lf %lf",&rx,&ry,&rz,&tau);
  21. printf("radii are %g %g %g %g\n",rx,ry,rz,tau);
  22. for(ib=0;ib<NB;ib++){
  23. sprintf(filename,"thermal_gauss%02d.dat",ib+1);
  24. fptr=fopen(filename,"w");
  25. for(iblankline=0;iblankline<3;iblankline++){
  26. fprintf(fptr,"line number %d, blah blah blah\n",iblankline);
  27. }
  28. fprintf(fptr,"%d %d %lf %lf\n",1,10000,0.0,0.0);
  29. for(i=0;i<10000;i++){
  30. id=(int)floor(NTYPES*ran2());
  31. if(mass[id]/temp>8.0){
  32. pmag2=0.0;
  33. gauss(&p[1],&p[2]);gauss(&p[3],&dummy);
  34. for(alpha=1;alpha<4;alpha++) {
  35. p[alpha]=sqrt(temp*mass[id])*p[alpha];
  36. pmag2=pmag2+p[alpha]*p[alpha];
  37. }
  38. p[0]=sqrt(pmag2+mass[id]*mass[id]);
  39. }
  40. else{
  41. do{
  42. pmag=-temp*log(ran2()*ran2()*ran2());
  43. p[0]=sqrt(mass[id]*mass[id]+pmag*pmag);
  44. weight=exp(-(p[0]-pmag)/temp);
  45. }while(ran2()>weight);
  46. //printf("id=%d pmag=%g\n",ident[id],pmag);
  47. cthet=1.0-2.0*ran2();
  48. sthet=sqrt(1.0-cthet*cthet);
  49. phi=2.0*pi*ran2();
  50. p[3]=pmag*cthet;
  51. p[1]=pmag*sthet*cos(phi);
  52. p[2]=pmag*sthet*sin(phi);
  53. }
  54. gauss(&r[0],&r[1]);gauss(&r[2],&r[3]);
  55. r[0]=tau*r[0];
  56. r[1]=rx*r[1];
  57. r[2]=ry*r[2];
  58. r[3]=rz*r[3];
  59. fprintf(fptr,"%d %d %9g %9g %9g %9g %9g %9g %9g %9g %9g\n",
  60. i+1,ident[id],p[1]/1000.0,p[2]/1000.0,p[3]/1000.0,p[0]/1000.0,
  61. mass[id]/1000.0,r[1],r[2],r[3],r[0]);
  62. }
  63. fclose(fptr);
  64. }
  65. return 0;
  66. }
  67. #define IM1 2147483563
  68. #define IM2 2147483399
  69. #define AM (1.0/IM1)
  70. #define IMM1 (IM1-1)
  71. #define IA1 40014
  72. #define IA2 40692
  73. #define IQ1 53668
  74. #define IQ2 52774
  75. #define IR1 12211
  76. #define IR2 3791
  77. #define NTAB 32
  78. #define NDIV (1+IMM1/NTAB)
  79. #define EPS 1.2e-7
  80. #define RNMX (1.0-EPS)
  81. double ran2() {
  82. long J;
  83. long K;
  84. static long IDUM2=123456789;
  85. static long IY=0;
  86. static long IV[NTAB];
  87. double randy;
  88. if(IDUM <=0){
  89. if(-(IDUM)<1) IDUM=1;
  90. else IDUM = -(IDUM);
  91. IDUM2=(IDUM);
  92. for (J=NTAB+7;J>=0;J--) {
  93. K=(IDUM)/IQ1;
  94. IDUM=IA1*(IDUM-K*IQ1)-K*IR1;
  95. if(IDUM<0) IDUM +=IM1;
  96. if(J<NTAB) IV[J]=IDUM;
  97. }
  98. IY=IV[0];
  99. //printf("Initialization Completed %d\n",IDUM);
  100. }
  101. K=(IDUM)/IQ1;
  102. IDUM=IA1*(IDUM-K*IQ1)-K*IR1;
  103. if(IDUM<0) IDUM+=IM1;
  104. K=IDUM2/IQ2;
  105. IDUM2=IA2*(IDUM2-K*IQ2)-K*IR2;
  106. if(IDUM2<0) IDUM2+=IM2;
  107. J=IY/NDIV;
  108. IY=IV[J]-IDUM2;
  109. IV[J]=IDUM;
  110. if(IY<1) IY+=IMM1;
  111. randy=AM*IY;
  112. if(randy>RNMX) return RNMX;
  113. else return randy;
  114. }
  115. //**********************************************************
  116. void gauss(double *g1,double *g2){
  117. double x,y,r2,r;
  118. do{
  119. x=1.0-2.0*ran2();
  120. y=1.0-2.0*ran2();
  121. r2=x*x+y*y;
  122. }while (r2>1.0);
  123. r=sqrt(r2);
  124. *g1=(x/r)*sqrt(-2.0*log(r2));
  125. *g2=(y/x)**g1;
  126. }