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 for(jj=0;jj<NTRAJ;jj++) 108 { 109 eng=engt; 110 111 engo = (eng+511.0)/(eng+1024.0); 112 sigma = p4*4.0e-18*z165*engo*engo/ 113 (eng+p3*1.0e-2*z133*pow(eng,0.5)); 114 lamda = lamda_a/sigma; 115 alpha = 7.0e-3/eng; 116 /*printf("engo,sigma,lamda,alpha %e %e %e %e \n", 117 engo,sigma,lamda,alpha); 118 */ 119 120 result = rand_mr(); 121 ran1=result/base; 122 step = -lamda*log(ran1); 123 /*printf("ran1,step %e %e\n",ran1,step); 124 */ 125 126 zcor=step; 127 cz=1; 128 129 de=step*p_sa*log(j_a*eng+1.0)/eng; 130 eng=eng+de; 131 /*printf("de,eng %e %e\n",de,eng); 132 */ 133 134 /*-------------------------------- 135 * now the loop for tracing the trajectory until it backscattered 136 * or it has energy too small to backscatter 137 */ 138 filtf=zcor>0 && eng>E_MIN && eng>(-de*zcor/step); 139 icic=filtf; 140 /* printf("\n"); 141 printf("icic,zcor,eng,de,step %d %e %e %e %e\n", 142 icic,zcor,eng,de,step); 143 */ 144 while(icic) 145 { 146 147 r_over_i=p1*100.0*pow(eng,zex)/zz+zz*zz*zz/eng/p2/1.0e5; 148 149 engo=(eng+511)/(eng+1024); 150 sigma = p4*4.0e-18*z165*engo*engo/ 151 (eng+p3*1.0e-2*z133*pow(eng,0.5)); 152 lamda= lamda_a/sigma; 153 alpha=7.0e-3/eng; 154 155 result = rand_mr(); 156 ran1=result/base; 157 result = rand_mr(); 158 ran2=result/base; 159 result = rand_mr(); 160 ran3=result/base; 161 result = rand_mr(); 162 ran4=result/base; 163 filtr=ran4<r_over_i/(r_over_i+1.0); 164 165 step = -lamda*log(ran1); 166 c_phi=filtr*(-2*alpha*ran2/(alpha-ran2+1)+1)+(1-filtr)*(1-2*ran2); 167 s_phi=sqrt(1.0-c_phi*c_phi); 168 psi=2*PI*ran3; 169 sz=sqrt(1.0-cz*cz); 170 cc = -s_phi*cos(psi)*sz+c_phi*cz; 171 172 zcor=zcor+step*cc; 173 de=step*p_sa*log(j_a*eng+1.0)/eng; 174 eng=eng+de; 175 cz=cc; 176 177 filtf=zcor>0 && eng>E_MIN && eng>(-de*zcor/step); 178 icic=filtf; 179 180 } /* end of the tracing process for this trajectory */ 181 182 if(zcor<0) iback++; 183 184 } /* end of the NTRAJ trajectories */ 185 186 #ifdef _DOPRINT_ 187 printf("point: engt %f, zz %f, A %f, rho %f\n",engt,zz,aa,rho); 188 printf("backscattering coeff for this point is %f \n",*backsc); 189 #endif 190 return (double)iback/(double)NTRAJ; 191 192 } /* end of back subroutine */ 193 /**************************************************************/ 194 195 196 /**********************************************************************/ 197 198 /**********************************************************************/ 199 /**********************************************************************/ 200 /********************** main function ***********************/ 201 /********************** ***********************/ 202 203 main(int argc, char ** argv) 204 { 205 int i,j,k,m; 206 int ii,jj,icic,icic2,icc; 207 double coe[NPAR]; 208 const char *fitfile, *montefile, *paramfile; 209 FILE *fpot, *fcoe, *fbac; 210 double point[POINT][4]; 211 double backsct[POINT],baexp[POINT],ajunk,bjunk; 212 double pointt[4],coet[NPAR],backcoe; 213 double rms,oldrms; 214 215 /* ---------------------------------------------------------- 216 * Read in the elements that need to be considered and the energy that 217 * need to be calculated for the back-scattering coefficient 218 */ 219 220 if (argc < 4) { 221 fprintf(stderr, "Not all files specified.\n Correct command line: %s fit.inp monte.dat param.dat\n", argv[0]); 222 exit(1); 223 } 224 225 fitfile = argv[1]; 226 montefile = argv[2]; 227 paramfile = argv[3]; 228 229 230 if ((fpot=fopen(fitfile,"r")) == NULL) 231 { 232 printf ("Error in opening %s\n", fitfile); 233 exit(1); 234 } 235 236 if ((fcoe=fopen(paramfile,"r")) == NULL) 237 { 238 printf ("Error in opening %s\n", paramfile); 239 exit(1); 240 } 241 242 if ((fbac=fopen(montefile,"r")) == NULL) 243 { 244 printf ("Error in opening %s\n", montefile); 245 exit(1); 246 } 247 /*-----------------------------------------------*/ 248 249 /*-----------------------------------------------*/ 250 /* Now, read in the data */ 251 252 for(i=0;i<POINT;i++) 253 { 254 for(j=0;j<4;j++) fscanf(fpot,"%lf",&point[i][j]); 255 for(j=0;j<4;j++) printf("%e ",point[i][j]); 256 printf("\n"); 257 } 258 printf("\n"); 259 260 for(i=0;i<NPAR;i++) 261 { 262 fscanf(fcoe,"%lf",&coe[i]); 263 printf("%e \n",coe[i]); 264 } 265 printf("\n"); 266 267 for(i=0;i<POINT;i++) 268 { 269 fscanf(fbac,"%lf%lf%lf",&ajunk,&bjunk,&baexp[i]); 270 printf("%e %e %e\n",ajunk,bjunk,baexp[i]); 271 } 272 printf("\n\n"); 273 fflush(stdout); 274 275 fclose(fpot); 276 fclose(fcoe); 277 fclose(fbac); 278 279 /****************************************************/ 280 /*----------------------- 281 * now calculate the POINT points of backscattering coefficients 282 */ 283 284 printf("first plot\n"); 285 printf("parameter set: \n"); 286 for(i=0;i<NPAR;i++) printf(" %e ",coe[i]); 287 printf("\n\n"); 288 fflush(stdout); 289 290 for(i=0;i<NPAR;i++) coet[i]=coe[i]; 291 292 #pragma omp parallel for private (i, pointt, backcoe) schedule(static) 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 #pragma omp parallel for private (i, pointt, backcoe) schedule(dynamic) 325 for(ii=0;ii<POINT;ii++) 326 { 327 for(i=0;i<4;i++) pointt[i]=point[ii][i]; 328 backcoe = back(pointt,coet); 329 backsct[ii]=backcoe; 330 fflush(stdout); 331 } 332 333 rms=rmse(backsct,baexp); 334 icic2=rms<oldrms; 335 336 if(icic2) 337 while(icic2) 338 { 339 icc=1; 340 oldrms=rms; 341 for(i=0;i<NPAR;i++) coe[i]=coet[i]; 342 printf("new parameter set: \n"); 343 for(i=0;i<NPAR;i++) printf(" %e ",coe[i]); 344 printf("\n"); 345 printf("corresponding rms: %f \n\n",oldrms); 346 fflush(stdout); 347 348 printf("increasing # %d \n",jj); 349 fflush(stdout); 350 351 for(i=0;i<NPAR;i++) coet[i]=coe[i]; 352 coet[jj]=coe[jj]*1.2; 353 354 #pragma omp parallel for private (i, pointt, backcoe) schedule(dynamic) 355 for(ii=0;ii<POINT;ii++) 356 { 357 for(i=0;i<4;i++) pointt[i]=point[ii][i]; 358 backcoe = back(pointt,coet); 359 backsct[ii]=backcoe; 360 fflush(stdout); 361 } 362 363 rms=rmse(backsct,baexp); 364 icic2=rms<oldrms; 365 } /* end of while and if */ 366 367 else 368 { 369 for(i=0;i<NPAR;i++) coet[i]=coe[i]; 370 coet[jj]=coe[jj]*0.8; 371 372 #pragma omp parallel for private (i, pointt, backcoe) schedule(dynamic) 373 for(ii=0;ii<POINT;ii++) 374 { 375 for(i=0;i<4;i++) pointt[i]=point[ii][i]; 376 backcoe = back(pointt,coet); 377 backsct[ii]=backcoe; 378 fflush(stdout); 379 } 380 381 rms=rmse(backsct,baexp); 382 icic2=rms<oldrms; 383 384 while(icic2) 385 { 386 icc=1; 387 oldrms=rms; 388 for(i=0;i<NPAR;i++) coe[i]=coet[i]; 389 printf("new parameter set: \n"); 390 for(i=0;i<NPAR;i++) printf(" %e ",coe[i]); 391 printf("\n"); 392 printf("corresponding rms: %f \n\n",oldrms); 393 fflush(stdout); 394 395 printf("decreasing # %d \n",jj); 396 fflush(stdout); 397 398 for(i=0;i<NPAR;i++) coet[i]=coe[i]; 399 coet[jj]=coe[jj]*0.8; 400 401 #pragma omp parallel for private (i, pointt, backcoe) schedule(dynamic) 402 for(ii=0;ii<POINT;ii++) 403 { 404 for(i=0;i<4;i++) pointt[i]=point[ii][i]; 405 backcoe = back(pointt,coet); 406 backsct[ii]=backcoe; 407 fflush(stdout); 408 } 409 410 rms=rmse(backsct,baexp); 411 icic2=rms<oldrms; 412 } /* end of while */ 413 } /* end of else */ 414 } /* end of for(jj=0... */ 415 416 icic=icc*(oldrms>ERRCRI); 417 418 } /* end of while(icic) */ 419 420 printf("\n\n***********************************************\n\n"); 421 printf("End of fitting process\n\n"); 422 printf("final parameter set: \n"); 423 for(i=0;i<NPAR;i++) printf(" %e ",coe[i]); 424 printf("\n\n Corresponding rmse: %f \n\n",oldrms); 425 printf("----------------------------------------\n\n"); 426 fflush(stdout); 427 428 }