1 #include <omp.h> 2 /* 3 * File: fit.c 4 * This is a general non-linear fitting program. It's used to a an 5 * optimal set of coefficients for the r_over _i ratio. 6 */ 7 8 #include <stdio.h> 9 #include <stdlib.h> 10 #include <math.h> 11 #include <sys/time.h> 12 #include "myrandom.h" 13 14 15 //#define _DOPRINT_ // #define NTRAJ 2500 //#define POINT 51 16 17 #define NTRAJ 750 18 #define E_MIN 0.05 19 #define POINT 51 20 #define NPAR 4 21 #define ERRCRI 0.01 22 #define PI 3.14159265358979 23 24 unsigned int INITIAL_SEED = 123; 25 26 /*---------------*/ 27 /* function 'rmse'. 28 * This function calculates the RMS error for the inputting 29 * bat and ba 30 * 'bat' is the back scattering coefficient of our own MC simulation 31 * 'ba' is the full quantum MC simulation and is used as experiment 32 * to which we are fitting to 33 */ 34 double rmse(double bat[], double ba[]); 35 36 /* function: rmse */ 37 double rmse(double bat[], double ba[]) 38 { 39 double rms; 40 int i; 41 42 rms=0.0; 43 for(i=0;i<POINT;i++) 44 { 45 rms += (bat[i]-ba[i])*(bat[i]-ba[i])/bat[i]/(1.0-bat[i]); 46 } 47 48 return (rms); 49 } 50 51 52 /*---------------*/ 53 /* private function 'back'. 54 * This function simulates the backscattering coefficient for the 55 * given element and energy 56 * array point stores the energy and element information in the 57 * sequence E,Z,A,rho (energy, atomic#, atomic weight, density) 58 * array coef is the array of the parameters (the parameters we are 59 * adjusting to make a good fit) 60 * backsc is the address of the output back scattering coefficient 61 */ 62 63 64 /* Function: backsc */ 65 double back(double point[], double coef[]) 66 { 67 int i,j,k,jj; 68 int icic,iback; 69 double eng; 70 double engo,sigma,lamda,alpha,r_over_i; 71 double ran1,ran2,ran3,ran4,step,zcor,cz,sz,de,c_phi,s_phi,psi,cc; 72 int filtr,filtf; 73 double base=(double)2147483647; 74 75 /*-----------------------*/ 76 /* Here is the section for setting up the some initial values for 77 * calculating sigma, alpha..., so that these values won't need to be 78 * recalculated for each trajectory. This section need modification if 79 * wemodify the formulas of alpha, r_over_i.. 80 */ 81 double p1=coef[0],p2=coef[1],p3=coef[2],p4=coef[3]; 82 double engt=point[0],zz=point[1],aa=point[2],rho=point[3]; 83 double z133=pow(zz,1.33),z165=pow(zz,1.65),zex=(1.0-zz/1000.0); 84 double lamda_a = aa/6.022e23/rho; 85 double j_j=0.001*(9.76*zz+58.5/pow(zz,0.19)); 86 double j_a=1.166/j_j; 87 double p_sa = -78500.0*rho*zz/aa; 88 /*-----------------------*/ 89 90 /*printf("p1,p2,p3,p4,p5,p6 %e %e %e %e %e %e\n",p1,p2,p3,p4,p5,p6); 91 printf("z133,z15 %e %e\n",z133,z15); 92 printf("lamda_a,j_j,j_a,p_sa, %e %e %e %e\n",lamda_a,j_j,j_a,p_sa); 93 */ 94 95 long int result; 96 97 srand_mr(INITIAL_SEED); 98 99 /*---------------------------------------------------------------- 100 * 101 * Now, we begin to calculate the back-scattering coefficient for 102 * this set of coefficient 103 */ 104 105 /* now, the loop for all NTRAJ trajectories */ 106 iback=0; // accumulator... 107 #pragma omp parallel for private (lamda, r_over_i, engo, cc, cz, result, eng, filtf, c_phi, filtr, psi, icic, s_phi, de, ran4, step, ran1, ran3, ran2, alpha, sz, zcor, sigma) schedule(static) 108 for(jj=0;jj<NTRAJ;jj++) 109 { 110 eng=engt; 111 112 engo = (eng+511.0)/(eng+1024.0); 113 sigma = p4*4.0e-18*z165*engo*engo/ 114 (eng+p3*1.0e-2*z133*pow(eng,0.5)); 115 lamda = lamda_a/sigma; 116 alpha = 7.0e-3/eng; 117 /*printf("engo,sigma,lamda,alpha %e %e %e %e \n", 118 engo,sigma,lamda,alpha); 119 */ 120 121 result = rand_mr(); 122 ran1=result/base; 123 step = -lamda*log(ran1); 124 /*printf("ran1,step %e %e\n",ran1,step); 125 */ 126 127 zcor=step; 128 cz=1; 129 130 de=step*p_sa*log(j_a*eng+1.0)/eng; 131 eng=eng+de; 132 /*printf("de,eng %e %e\n",de,eng); 133 */ 134 135 /*-------------------------------- 136 * now the loop for tracing the trajectory until it backscattered 137 * or it has energy too small to backscatter 138 */ 139 filtf=zcor>0 && eng>E_MIN && eng>(-de*zcor/step); 140 icic=filtf; 141 /* printf("\n"); 142 printf("icic,zcor,eng,de,step %d %e %e %e %e\n", 143 icic,zcor,eng,de,step); 144 */ 145 while(icic) 146 { 147 148 r_over_i=p1*100.0*pow(eng,zex)/zz+zz*zz*zz/eng/p2/1.0e5; 149 150 engo=(eng+511)/(eng+1024); 151 sigma = p4*4.0e-18*z165*engo*engo/ 152 (eng+p3*1.0e-2*z133*pow(eng,0.5)); 153 lamda= lamda_a/sigma; 154 alpha=7.0e-3/eng; 155 156 result = rand_mr(); 157 ran1=result/base; 158 result = rand_mr(); 159 ran2=result/base; 160 result = rand_mr(); 161 ran3=result/base; 162 result = rand_mr(); 163 ran4=result/base; 164 filtr=ran4<r_over_i/(r_over_i+1.0); 165 166 step = -lamda*log(ran1); 167 c_phi=filtr*(-2*alpha*ran2/(alpha-ran2+1)+1)+(1-filtr)*(1-2*ran2); 168 s_phi=sqrt(1.0-c_phi*c_phi); 169 psi=2*PI*ran3; 170 sz=sqrt(1.0-cz*cz); 171 cc = -s_phi*cos(psi)*sz+c_phi*cz; 172 173 zcor=zcor+step*cc; 174 de=step*p_sa*log(j_a*eng+1.0)/eng; 175 eng=eng+de; 176 cz=cc; 177 178 filtf=zcor>0 && eng>E_MIN && eng>(-de*zcor/step); 179 icic=filtf; 180 181 } /* end of the tracing process for this trajectory */ 182 183 if(zcor<0) iback++; 184 185 } /* end of the NTRAJ trajectories */ 186 187 #ifdef _DOPRINT_ 188 printf("point: engt %f, zz %f, A %f, rho %f\n",engt,zz,aa,rho); 189 printf("backscattering coeff for this point is %f \n",*backsc); 190 #endif 191 return (double)iback/(double)NTRAJ; 192 193 } /* end of back subroutine */ 194 /**************************************************************/ 195 196 197 /**********************************************************************/ 198 199 /**********************************************************************/ 200 /**********************************************************************/ 201 /********************** main function ***********************/ 202 /********************** ***********************/ 203 204 main(int argc, char ** argv) 205 { 206 int i,j,k,m; 207 int ii,jj,icic,icic2,icc; 208 double coe[NPAR]; 209 const char *fitfile, *montefile, *paramfile; 210 FILE *fpot, *fcoe, *fbac; 211 double point[POINT][4]; 212 double backsct[POINT],baexp[POINT],ajunk,bjunk; 213 double pointt[4],coet[NPAR],backcoe; 214 double rms,oldrms; 215 216 /* ---------------------------------------------------------- 217 * Read in the elements that need to be considered and the energy that 218 * need to be calculated for the back-scattering coefficient 219 */ 220 221 if (argc < 4) { 222 fprintf(stderr, "Not all files specified.\n Correct command line: %s fit.inp monte.dat param.dat\n", argv[0]); 223 exit(1); 224 } 225 226 fitfile = argv[1]; 227 montefile = argv[2]; 228 paramfile = argv[3]; 229 230 231 if ((fpot=fopen(fitfile,"r")) == NULL) 232 { 233 printf ("Error in opening %s\n", fitfile); 234 exit(1); 235 } 236 237 if ((fcoe=fopen(paramfile,"r")) == NULL) 238 { 239 printf ("Error in opening %s\n", paramfile); 240 exit(1); 241 } 242 243 if ((fbac=fopen(montefile,"r")) == NULL) 244 { 245 printf ("Error in opening %s\n", montefile); 246 exit(1); 247 } 248 /*-----------------------------------------------*/ 249 250 /*-----------------------------------------------*/ 251 /* Now, read in the data */ 252 253 for(i=0;i<POINT;i++) 254 { 255 for(j=0;j<4;j++) fscanf(fpot,"%lf",&point[i][j]); 256 for(j=0;j<4;j++) printf("%e ",point[i][j]); 257 printf("\n"); 258 } 259 printf("\n"); 260 261 for(i=0;i<NPAR;i++) 262 { 263 fscanf(fcoe,"%lf",&coe[i]); 264 printf("%e \n",coe[i]); 265 } 266 printf("\n"); 267 268 for(i=0;i<POINT;i++) 269 { 270 fscanf(fbac,"%lf%lf%lf",&ajunk,&bjunk,&baexp[i]); 271 printf("%e %e %e\n",ajunk,bjunk,baexp[i]); 272 } 273 printf("\n\n"); 274 fflush(stdout); 275 276 fclose(fpot); 277 fclose(fcoe); 278 fclose(fbac); 279 280 /****************************************************/ 281 /*----------------------- 282 * now calculate the POINT points of backscattering coefficients 283 */ 284 285 printf("first plot\n"); 286 printf("parameter set: \n"); 287 for(i=0;i<NPAR;i++) printf(" %e ",coe[i]); 288 printf("\n\n"); 289 fflush(stdout); 290 291 for(i=0;i<NPAR;i++) coet[i]=coe[i]; 292 293 for(ii=0;ii<POINT;ii++) 294 { 295 for(i=0;i<4;i++) pointt[i]=point[ii][i]; 296 backcoe = back(pointt,coet); 297 backsct[ii]=backcoe; 298 fflush(stdout); 299 } 300 301 rms=rmse(backsct,baexp); 302 oldrms=rms; 303 icic=oldrms>ERRCRI; 304 printf("icic %d, rms %f \n\n",icic,oldrms); 305 fflush(stdout); 306 307 /*---------------------------------*/ 308 /* Now we begin the loop for automated fitting */ 309 310 while(icic) 311 { 312 313 /*----------------------------------------------------*/ 314 icc=0; 315 316 for(jj=0;jj<NPAR;jj++) 317 { 318 printf("Changing # %d \n\n",jj); 319 fflush(stdout); 320 321 for(i=0;i<NPAR;i++) coet[i]=coe[i]; 322 coet[jj]=coe[jj]*1.2; 323 324 for(ii=0;ii<POINT;ii++) 325 { 326 for(i=0;i<4;i++) pointt[i]=point[ii][i]; 327 backcoe = back(pointt,coet); 328 backsct[ii]=backcoe; 329 fflush(stdout); 330 } 331 332 rms=rmse(backsct,baexp); 333 icic2=rms<oldrms; 334 335 if(icic2) 336 while(icic2) 337 { 338 icc=1; 339 oldrms=rms; 340 for(i=0;i<NPAR;i++) coe[i]=coet[i]; 341 printf("new parameter set: \n"); 342 for(i=0;i<NPAR;i++) printf(" %e ",coe[i]); 343 printf("\n"); 344 printf("corresponding rms: %f \n\n",oldrms); 345 fflush(stdout); 346 347 printf("increasing # %d \n",jj); 348 fflush(stdout); 349 350 for(i=0;i<NPAR;i++) coet[i]=coe[i]; 351 coet[jj]=coe[jj]*1.2; 352 353 for(ii=0;ii<POINT;ii++) 354 { 355 for(i=0;i<4;i++) pointt[i]=point[ii][i]; 356 backcoe = back(pointt,coet); 357 backsct[ii]=backcoe; 358 fflush(stdout); 359 } 360 361 rms=rmse(backsct,baexp); 362 icic2=rms<oldrms; 363 } /* end of while and if */ 364 365 else 366 { 367 for(i=0;i<NPAR;i++) coet[i]=coe[i]; 368 coet[jj]=coe[jj]*0.8; 369 370 for(ii=0;ii<POINT;ii++) 371 { 372 for(i=0;i<4;i++) pointt[i]=point[ii][i]; 373 backcoe = back(pointt,coet); 374 backsct[ii]=backcoe; 375 fflush(stdout); 376 } 377 378 rms=rmse(backsct,baexp); 379 icic2=rms<oldrms; 380 381 while(icic2) 382 { 383 icc=1; 384 oldrms=rms; 385 for(i=0;i<NPAR;i++) coe[i]=coet[i]; 386 printf("new parameter set: \n"); 387 for(i=0;i<NPAR;i++) printf(" %e ",coe[i]); 388 printf("\n"); 389 printf("corresponding rms: %f \n\n",oldrms); 390 fflush(stdout); 391 392 printf("decreasing # %d \n",jj); 393 fflush(stdout); 394 395 for(i=0;i<NPAR;i++) coet[i]=coe[i]; 396 coet[jj]=coe[jj]*0.8; 397 398 for(ii=0;ii<POINT;ii++) 399 { 400 for(i=0;i<4;i++) pointt[i]=point[ii][i]; 401 backcoe = back(pointt,coet); 402 backsct[ii]=backcoe; 403 fflush(stdout); 404 } 405 406 rms=rmse(backsct,baexp); 407 icic2=rms<oldrms; 408 } /* end of while */ 409 } /* end of else */ 410 } /* end of for(jj=0... */ 411 412 icic=icc*(oldrms>ERRCRI); 413 414 } /* end of while(icic) */ 415 416 printf("\n\n***********************************************\n\n"); 417 printf("End of fitting process\n\n"); 418 printf("final parameter set: \n"); 419 for(i=0;i<NPAR;i++) printf(" %e ",coe[i]); 420 printf("\n\n Corresponding rmse: %f \n\n",oldrms); 421 printf("----------------------------------------\n\n"); 422 fflush(stdout); 423 424 }