1/******************************************************************************* 2 * 3 * McStas, neutron ray-tracing package 4 * Copyright (C) 1997-2008, All rights reserved 5 * Risoe National Laboratory, Roskilde, Denmark 6 * Institut Laue Langevin, Grenoble, France 7 * 8 * Component: ViewModISIS 9 * 10 * 11 * %I 12 * 13 * Modification of ViewModerator4 component written by S. Ansell in 2015 14 * 15 * G. Skoro 16 * Date: February 2016 17 * Version: 1st Draft 18 * Origin: ISIS 19 * Tested with McStas 2.0 and 2.2 (Windows) 20 * 21 * 22 * %D 23 * Produces a neutron distribution at the ISIS TS1 or TS2 beamline shutter front face (or corresponding moderator) position. 24 * The Face argument determines which TS1 or TS2 beamline is to be sampled by using corresponding Etable file. 25 * Neutrons are created having a range of energies determined by the E0 and E1 arguments. 26 * Trajectories are produced such that they pass through the shutter front face (RECOMENDED) or moderator face (defined by 27 * xw and yh) and a focusing rectangle (defined by focus_xw, focus_yh and dist). 28 * 29 * 30 * %P 31 * INPUT PARAMETERS: 32 * 33 * Face: [string] Etable filename 34 * E0: [meV] Lower edge of energy distribution 35 * E1: [meV] Upper edge of energy distribution 36 * modPosition: [int] Defines the initial surface for neutron distribution production. Possible values = 0 (at moderator) or 1 (at shutter insert) 37 * xw: [m] Moderator width 38 * yh: [m] Moderator height 39 * focus_xw: [m] Width of focusing rectangle 40 * focus_yh: [m] Height of focusing rectangle 41 * dist: [m] Distance from source surface to the focusing rectangle 42 * verbose: [int] Flag to output debugging information 43 * beamcurrent [uA] ISIS beam current 44 * 45 * Example 1: ViewModISISver1(Face="TS1_S04_Merlin.mcstas", E0 = E_min, E1 = E_max, 46 * dist = 1.7, focus_xw = 0.094, focus_yh = 0.094, modPosition=0, 47 * xw=0.12,yh=0.12) 48 * 49 * MERLIN simulation. 50 * modPosition=0 so initial surface is at (near to) the moderator face. 51 * In this example, focus_xw and focus_yh are chosen to be identical to the shutter opening dimension. 52 * dist = 1.7 is the 'real' distance to the shutter front face => This is TimeOffset value (=170 [cm]) 53 * from TS1_S04_Merlin.mcstas file. 54 * 55 * 56 * Example 2: ViewModISISver1(Face="Let_timeTest_155.mcstas", E0 = E_min, E1 = E_max, 57 * modPosition=1, xw=0.196, yh = 0.12, focus_xw = 0.04, focus_yh = 0.094, dist = 0.5) 58 * 59 * LET simulation. 60 * modPosition=1 so initial surface is at front face of shutter insert. 61 * 62 * (IMPORTANT) If modPosition=1, the xw and yh values are obsolete. The dimensions of 63 * initial surface are automatically calculated using "RDUM" values in Let_timeTest_155.mcstas file. 64 * 65 * In this example, focus_xw and focus_yh are arbitrary chosen to be identical to the shutter opening dimension. 66 * 67 * 68 * 69 * N.B. Absolute normalization: The result of the Mc-Stas simulation will show neutron intensity for beam current of 1 uA. 70 * 71 * 72 * %E 73 *******************************************************************************/ 74 75DEFINE COMPONENT ViewModISIS 76SETTING PARAMETERS (string Face="TS1_S04_Merlin.mcstas",E0, E1,modPosition=0,xw=0.12,yh=0.12,focus_xw=0.094,focus_yh=0.094,dist=1.7,int verbose=0, beamcurrent=1) 77DECLARE 78%{ 79 80#include <ctype.h> 81 82typedef struct 83{ 84 int nEnergy; ///< Number of energy bins 85 int nTime; ///< number of time bins 86 87 double XAxis; 88 double ZAxis; 89 90 double rdumMid; ///< tally time Window mid point 91 double timeOffset; ///< Time separation 92 double* TimeBin; ///< Time bins 93 double* EnergyBin; ///< Energy bins 94 95 double** Flux; ///< Flux per bin (integrated) 96 double* EInt; ///< Integrated Energy point 97 double Total; ///< Integrated Total 98 99} Source; 100 101 /* New functions */ 102 103 int cmdnumberD(char *,double*); 104 int cmdnumberI(char *,int*,const int); 105 double polInterp(double*,double*,int,double); 106 int findRDUM(char*); 107 int findTimeOffset(char*,double*); 108 FILE* openFile(char*); 109 FILE* openFileTest(char*); 110 void processHeader(FILE*); 111 void readHtable(FILE*,const double,const double); 112 int timeStart(char*); 113 int timeEnd(char*); 114 int energyBin(char*,double,double,double*,double*); 115 int notComment(char*); 116 void orderEnergy(double*,double*); 117 void cutToNumber(char*); 118 double convertEnergy(double); 119 double EtoLambda(double); 120 double calcRDum(double*); 121 122 double strArea(); 123 124 /* global variables */ 125 126 double p_in; /* Polorization term (from McSTAS) */ 127 int Tnpts; /* Number of points in parameteriation */ 128 129 double scaleSize; /* correction for the actual area of the moderator viewed */ 130 131 double angleArea; /* Area seen by the window */ 132 133 double Nsim; /* Total number of neutrons to be simulated */ 134 135 int Ncount; /* Number of neutron simulate so far*/ 136 Source TS; 137 138 /* runtime variables*/ 139 double fullAngle; 140 double rtE0,rtE1; /* runtime Energy minima and maxima so we can use angstroms as negative input */ 141 142 143double** 144matrix(const int m,const int n) 145 /*! 146 Determine a double matrix 147 */ 148{ 149 int i; 150 double* pv; 151 double** pd; 152 153 if (m<1) return 0; 154 if (n<1) return 0; 155 pv = (double*) malloc(m*n*sizeof(double)); 156 pd = (double**) malloc(m*sizeof(double*)); 157 if (!pd) 158 { 159 fprintf(stderr,"No room for matrix!\n"); 160 exit(1); 161 } 162 for (i=0;i<m;i++) 163 pd[i]=pv + (i*n); 164 return pd; 165} 166 167 168double 169polInterp(double* X,double* Y,int Psize,double Aim) 170 /*! 171 returns the interpolated polynomial between Epnts 172 and the integration 173 \param X :: X coordinates 174 \param Y :: Y coordinates 175 \param Psize :: number of valid point in array to use 176 \param Aim :: Aim point to intepolate result (X coordinate) 177 \returns Energy point 178 */ 179{ 180 double out,errOut; /* out put variables */ 181 double C[Psize],D[Psize]; 182 double testDiff,diff; 183 184 double w,den,ho,hp; /* intermediate variables */ 185 int i,m,ns; 186 187 188 ns=0; 189 diff=fabs(Aim-X[0]); 190 C[0]=Y[0]; 191 D[0]=Y[0]; 192 for(i=1;i<Psize;i++) 193 { 194 testDiff=fabs(Aim-X[i]); 195 if (diff>testDiff) 196 { 197 ns=i; 198 diff=testDiff; 199 } 200 C[i]=Y[i]; 201 D[i]=Y[i]; 202 } 203 204 out=Y[ns]; 205 ns--; /* Now can be -1 !!!! */ 206 207 for(m=1;m<Psize;m++) 208 { 209 for(i=0;i<Psize-m;i++) 210 { 211 ho=X[i]-Aim; 212 hp=X[i+m]-Aim; 213 w=C[i+1]-D[i]; 214 /* den=ho-hp; -- test !=0.0 */ 215 den=w/(ho-hp); 216 D[i]=hp*den; 217 C[i]=ho*den; 218 } 219 220 errOut= (2*(ns+1)<(Psize-m)) ? C[ns+1] : D[ns--]; 221 out+=errOut; 222 } 223 return out; 224} 225 226int 227binSearch(int Npts,double* AR,double V) 228 /*! 229 Object is to find the point in 230 array AR, closest to the value V 231 Checked for ordered array returns lower of backeting objects 232 */ 233{ 234 int klo,khi,k; 235 if (Npts<=0) 236 return 0; 237 if (V>AR[Npts-1]) 238 return Npts; 239 240 if(AR[0]>0.0)AR[0]=0.0; 241 242 if (V<AR[0]) 243 { 244 // if(AR[0]>0.0)AR[0]=0.0; 245 fprintf(stderr,"here"); 246 return 0; 247 } 248 klo=0; 249 khi= Npts-1; 250 while (khi-klo >1) 251 { 252 k=(khi+klo) >> 1; // quick division by 2 253 if (AR[k]>V) 254 khi=k; 255 else 256 klo=k; 257 } 258 return khi; 259} 260 261int 262cmdnumberD(char *mc,double* num) 263 /*! 264 \returns 1 on success 0 on failure 265 */ 266{ 267 int i,j; 268 char* ss; 269 char **endptr; 270 double nmb; 271 int len; 272 273 len=strlen(mc); 274 j=0; 275 276 for(i=0;i<len && mc[i] && 277 (mc[i]=='\t' || mc[i]==' ' || mc[i]==',');i++); 278 if(i==len || !mc[i]) return 0; 279 ss=malloc(sizeof(char)*(len+1)); 280 281 for(;i<len && mc[i]!='\n' && mc[i] 282 && mc[i]!='\t' && mc[i]!=' ' && mc[i]!=',';i++) 283 { 284 ss[j]=mc[i]; 285 j++; 286 } 287 if (!j) 288 { 289 free(ss); 290 return 0; //This should be impossible 291 } 292 ss[j]=0; 293 endptr=malloc(sizeof(char*)); 294 nmb = strtod(ss,endptr); 295 if (*endptr != ss+j) 296 { 297 free(endptr); 298 free(ss); 299 return 0; 300 } 301 *num = (double) nmb; 302 for(j=0;j<i && mc[j];j++) 303 mc[j]=' '; 304 free(endptr); 305 free(ss); 306 return 1; 307} 308 309int 310notComment(char* Line) 311 /*! 312 \returns 0 on a comment, 1 on a non-comment 313 */ 314{ 315 int len,i; 316 317 len=strlen(Line); 318 for(i=0;i<len && isspace(Line[i]);i++); 319 320 if (!Line[i] || Line[i]=='c' || Line[i]=='C' || 321 Line[i]=='!' || Line[i]=='#') 322 return 0; 323 return 1; 324} 325 326void 327orderEnergy(double* A,double *B) 328{ 329 double tmp; 330 if (*A>*B) 331 { 332 tmp=*A; 333 *A=*B; 334 *B=tmp; 335 } 336 return; 337} 338 339int 340timeStart(char* Line) 341 /*! 342 Search for a word time at the start of 343 the line. 344 \param Line :: Line to search 345 \returns 1 on success 0 on failure 346 */ 347{ 348 int len,i; 349 350 len=strlen(Line); 351 for(i=0;i<len && isspace(Line[i]);i++); 352 if (len-i<4) return 0; 353 return (strncmp(Line+i,"time",4)) ? 0 : 1; 354} 355 356void 357cutToNumber(char* Line) 358{ 359 int len,i; 360 361 len=strlen(Line); 362 for(i=0;i<len && !isdigit(Line[i]) && Line[i]!='-';i++) 363 Line[i]=' '; 364 return; 365} 366 367int 368findTimeOffset(char* Line,double* TO) 369 /*! 370 Determine the time offset from the file 371 \return 1 on success 372 */ 373{ 374 int len,i; 375 double D; 376 377 len=strlen(Line); 378 for(i=0;i<len && Line[i]!='T';i++); 379 if (len-i<11) return 0; 380 381 if (strncmp(Line+i,"TimeOffset",11) && 382 cmdnumberD(Line+i+11,&D)) 383 { 384 *TO=D/100.0; 385 return 1; 386 } 387 return 0; 388} 389 390int 391findRDUM(char* Line) 392 /*! 393 Search for a word rdum at the start of 394 the line. 395 \param Line :: Line to search 396 \returns 1 on success 0 on failure 397 */ 398{ 399 int len,i; 400 401 len=strlen(Line); 402 for(i=0;i<len && Line[i]!='r';i++); 403 if (len-i<4) return 0; 404 return (strncmp(Line+i,"rdum",4)) ? 0 : 1; 405} 406 407double 408convertEnergy(double E) 409 /*! 410 Convert the energy from eV [not change] or 411 from negative -ve which is angstrom 412 */ 413{ 414 return (E>0.0) ? E : 81.793936/(E*E); 415} 416 417double 418EtoLambda(double E) 419 /*! 420 Convert the energy from eV [not change] o 421 to lambda [A] 422 */ 423{ 424 return sqrt(81.793936/E); 425} 426 427 428int 429timeEnd(char* Line) 430 /*! 431 Search for a word time at the start of 432 the line. 433 \param Line :: Line to search 434 \returns 1 on success 0 on failure 435 */ 436{ 437 int len,i; 438 439 len=strlen(Line); 440 for(i=0;i<len && isspace(Line[i]);i++); 441 if (len-i<5) return 0; 442 return (strncmp(Line+i,"total",5)) ? 0 : 1; 443} 444 445int 446energyBin(char* Line,double Einit,double Eend,double* Ea,double* Eb) 447 /*! 448 Search for a word "energy bin:" at the start of 449 the line. Then separte off the energy bin values 450 \param Line :: Line to search 451 \param Ea :: first energy bin [meV] 452 \param Eb :: second energy bin [meV] 453 \returns 1 on success 0 on failure 454 */ 455{ 456 int len,i; 457 double A,B; 458 459 len=strlen(Line); 460 for(i=0;i<len && isspace(Line[i]);i++); 461 if (len-i<11) return 0; 462 463 464 if (strncmp(Line+i,"energy bin:",11)) 465 return 0; 466 467 i+=11; 468 if (!cmdnumberD(Line+i,&A)) 469 return 0; 470 // remove 'to' 471 for(;i<len-1 && Line[i]!='o';i++); 472 i++; 473 if (!cmdnumberD(Line+i,&B)) 474 return 0; 475 A*=1e9; 476 B*=1e9; 477 *Ea=A; 478 *Eb=B; 479 if (*Eb>Einit && *Ea<Eend) 480 return 1; 481 return 0; 482} 483 484double 485calcFraction(double EI,double EE,double Ea,double Eb) 486 /*! 487 Calculate the fraction of the bin between Ea -> Eb 488 that is encompassed by EI->EE 489 */ 490{ 491 double frac; 492 double dRange; 493 494 if (EI>Eb) 495 return 0.0; 496 if (EE<Ea) 497 return 0.0; 498 499 dRange=Eb-Ea; 500 frac=(EI>Ea) ? (Eb-EI)/dRange : 1.0; 501 502 503 frac-=(EE<Eb) ? (Eb-EE)/dRange : 0.0; 504 505 // if(frac != 1.0) 506 // fprintf(stderr,"frac %g, Ea %g,Eb %g, EI %g, EE %g\n",frac,Ea,Eb,EI,EE); 507 508 return frac; 509} 510 511double 512calcRDum(double* RPts) 513 /*! 514 Caluclate the mean distance from the window centre point 515 */ 516{ 517 double sum,zMid; 518 int i; 519 extern Source TS; 520 521 sum=0.0; 522 for(i=0;i<4;i++) 523 { 524 fprintf(stderr,"RDUM %g %g %g\n",RPts[i*3],RPts[i*3+1],RPts[i*3+2]); 525 sum+=sqrt(RPts[i*3]*RPts[i*3]+ 526 RPts[i*3+1]*RPts[i*3+1]); 527 } 528 sum/=400.0; // Convert to metres from cm 529 return sum; 530} 531 532void 533processHeader(FILE* TFile) 534 /*! 535 Determine the window and the time offset 536 \param TFile :: file [rewound on exit] 537 */ 538{ 539 char ss[255]; /* BIG space for line */ 540 int rdumCnt; 541 double Ea,Eb; 542 double RPts[12]; 543 int timeOffsetFlag; 544 int i,j; 545 double D; 546 547 548 extern Source TS; 549 550 rdumCnt=0; 551 timeOffsetFlag=0; 552 553 while(fgets(ss,255,TFile)) 554 { 555 if (!timeOffsetFlag) 556 timeOffsetFlag=findTimeOffset(ss,&TS.timeOffset); 557 558 if (!rdumCnt) 559 rdumCnt=findRDUM(ss); 560 if (rdumCnt && rdumCnt<5) 561 { 562 cutToNumber(ss); 563 for(i=0;i<3 && cmdnumberD(ss,&D);i++) 564 RPts[(rdumCnt-1)*3+i]=D; 565 rdumCnt++; 566 } 567 // EXIT CONDITION: 568 if (rdumCnt*timeOffsetFlag==5) 569 { 570 for(j=0;j<3;j++) 571 { 572 TS.ZAxis+=(RPts[3+j]-RPts[j])*(RPts[3+j]-RPts[j]); 573 TS.XAxis+=(RPts[6+j]-RPts[3+j])*(RPts[6+j]-RPts[3+j]); 574 } 575 576 577 TS.ZAxis=sqrt(TS.ZAxis)/100.0; // convert to metres from cm 578 TS.XAxis=sqrt(TS.XAxis)/100.0; 579 if (!modPosition) 580 { 581 TS.ZAxis=yh; 582 TS.XAxis=xw; 583 } 584 fprintf(stderr,"Time off sec == %g \n",TS.timeOffset); 585 fprintf(stderr,"mod size z == %g \n",TS.ZAxis); 586 // TS.rdumMid=calcRDum(RPts); 587 TS.rdumMid=TS.timeOffset; // Goran 588 return; 589 } 590 } 591 return; 592} 593 594 595 596void 597readHtable(FILE* TFile,const double Einit,const double Eend) 598/*! 599 Process a general h.o file to create an integrated 600 table of results from Einit -> Eend 601 \param Einit :: inital Energy 602 \param Eend :: final energy 603*/ 604{ 605 char ss[255]; /* BIG space for line */ 606 double Ea,Eb; 607 double T,D; 608 double Efrac; // Fraction of an Energy Bin 609 int Ftime; // time Flag 610 int eIndex; // energy Index 611 int tIndex; // time Index 612 double Tsum; // Running integration 613 double Efraction; // Amount to use for an energy/time bin 614 615 extern Source TS; 616 617 int i; 618 /*! 619 Status Flag:: 620 Ftime=1 :: [time ] Reading Time : Data : Err [Exit on Total] 621 622 Double Read File to determine how many bins and 623 memery size 624 */ 625 Ea=0.0; 626 Eb=0.0; 627 fprintf(stderr,"Energy == %g %g\n",Einit,Eend); 628 eIndex= -1; 629 Ftime=0; 630 tIndex=0; 631 TS.nTime=0; 632 TS.nEnergy=0; 633 634 processHeader(TFile); 635 636 // Read file and get time bins: 637 while(fgets(ss,255,TFile) && Eend>Ea) 638 { 639 if (notComment(ss)) 640 { 641 if (!Ftime) 642 { 643 // find :: energy bin: <Number> to <Number> 644 if (energyBin(ss,Einit,Eend,&Ea,&Eb)) 645 { 646 if (eIndex==0) 647 TS.nTime=tIndex; 648 eIndex++; 649 } 650 else if (timeStart(ss)) // find :: <spc> time 651 { 652 Ftime=1; 653 tIndex=0; 654 } 655 } 656 else // In the time data section 657 { 658 if (timeEnd(ss)) // Found "total" 659 Ftime=0; 660 else 661 { 662 // Need to read the line in the case of first run 663 if (TS.nTime==0) 664 { 665 if (cmdnumberD(ss,&T) && 666 cmdnumberD(ss,&D)) 667 tIndex++; 668 } 669 } 670 } 671 } 672 } 673 // 674 // Plus 2 since we have a 0 counter and we have missed the last line. 675 // 676 TS.nEnergy=eIndex+2; 677 if (!TS.nTime && tIndex) // SMALL TEST if reading first energy bin 678 TS.nTime=tIndex; 679 // printf("tIndex %d %d %d %d \n",tIndex,eIndex,TS.nEnergy,TS.nTime); 680 681 /* SECOND TIME THROUGH:: */ 682 rewind(TFile); 683 684 TS.Flux=matrix(TS.nEnergy,TS.nTime); 685 TS.EInt=(double*) malloc(TS.nEnergy*sizeof(double)); // Runing integral of energy 686 TS.TimeBin=(double*) malloc(TS.nTime*sizeof(double)); 687 TS.EnergyBin=(double*) malloc(TS.nEnergy*sizeof(double)); 688 if (verbose) fprintf(stderr,"NUMBER %d %d \n",TS.nEnergy,TS.nTime); 689 690 Tsum=0.0; 691 Ea=0.0; 692 Eb=0.0; 693 eIndex=-1; 694 Ftime=0; 695 tIndex=0; 696 TS.EInt[0]=0.0; 697 698 // Read file and get time bins 699 while(fgets(ss,255,TFile) && Eend>Ea) 700 { 701 if (notComment(ss)) 702 { 703 if (!Ftime) 704 { 705 if (energyBin(ss,Einit,Eend,&Ea,&Eb)) 706 { 707 eIndex++; 708 TS.EnergyBin[eIndex]=(Einit>Ea) ? Einit : Ea; 709 Efraction=calcFraction(Einit,Eend,Ea,Eb); 710 Ftime++; 711 } 712 } 713 else if (Ftime==1) 714 { 715 if (timeStart(ss)) 716 { 717 Ftime=2; 718 tIndex=0; 719 } 720 } 721 else // In the time section 722 { 723 if (timeEnd(ss)) // Found "total" 724 { 725 Ftime=0; 726 TS.EInt[eIndex+1]=Tsum; 727 728 if (verbose) fprintf(stderr,"Energy %g %g \n",EtoLambda(Ea),(TS.EInt[eIndex+1]-TS.EInt[eIndex])*3.744905847e14); 729 } 730 else 731 { 732 // Need to read the line in the case of first run 733 if (cmdnumberD(ss,&T) && 734 cmdnumberD(ss,&D)) 735 { 736 TS.TimeBin[tIndex]=T/1e8; // convert Time into second (from shakes) 737 Tsum+=D*Efraction; 738 TS.Flux[eIndex][tIndex]=Tsum; 739 tIndex++; 740 } 741 } 742 } 743 } 744 } 745 TS.EnergyBin[eIndex+1]=Eend; 746 TS.Total=Tsum; 747 748 // printf("tIndex %d %d %d \n",tIndex,eIndex,TS.nTime); 749 // printf("Tsum %g \n",Tsum); 750 // fprintf(stderr,"ebin1 ebinN %g %g\n",TS.EnergyBin[0],TS.EnergyBin[TS.nEnergy-1]); 751 752 return; 753} 754 755void 756getPoint(double* TV,double* EV,double* lim1, double* lim2) 757 /*! 758 Calculate the Time and Energy 759 by sampling the file. 760 Uses TS table to find the point 761 \param TV :: 762 \param EV :: 763 \param lim1 :: 764 \param lim2 :: 765 */ 766{ 767 int i; 768 769 extern Source TS; 770 double R0,R1,R,Rend; 771 int Epnt; ///< Points to the next higher index of the neutron integral 772 int Tpnt; 773 int iStart,iEnd; 774 double TRange,Tspread; 775 double Espread,Estart; 776 double *EX; 777 778 // So that lowPoly+highPoly==maxPoly 779 const int maxPoly=6; 780 const int highPoly=maxPoly/2; 781 const int lowPoly=maxPoly-highPoly; 782 783 // static int testVar=0; 784 785 R0=rand01(); 786 Rend=R=TS.Total*R0; 787 // This gives Eint[Epnt-1] > R > Eint[Epnt] 788 Epnt=binSearch(TS.nEnergy-1,TS.EInt,R); 789 Tpnt=binSearch(TS.nTime-1,TS.Flux[Epnt-1],R); 790 791 if(Epnt<0 || R<TS.Flux[Epnt-1][Tpnt-1] || R >TS.Flux[Epnt-1][Tpnt] ) 792 { 793 fprintf(stderr,"outside bin limits Tpnt/Epnt problem %12.6e %12.6e %12.6e \n", 794 TS.Flux[Epnt-1][Tpnt-1],R,TS.Flux[Epnt-1][Tpnt]); 795 } 796 797 if(Epnt == 0) 798 { 799 Estart=0.0; 800 Espread=TS.EInt[0]; 801 *EV=TS.EnergyBin[1]; 802 } 803 else 804 { 805 Estart=TS.EInt[Epnt-1]; 806 Espread=TS.EInt[Epnt]-TS.EInt[Epnt-1]; 807 *EV=(Epnt>TS.nEnergy-1) ? TS.EnergyBin[Epnt+1] : TS.EnergyBin[Epnt]; 808 } 809 810 if (Tpnt==0 || Epnt==0 || Tpnt==TS.nTime) 811 { 812 fprintf(stderr,"BIG ERROR WITH Tpnt: %d and Epnt: %d\n",Tpnt,Epnt); 813 exit(1); 814 } 815 *TV=TS.TimeBin[Tpnt-1]; 816 TRange=TS.TimeBin[Tpnt]-TS.TimeBin[Tpnt-1]; 817 Tspread=TS.Flux[Epnt-1][Tpnt]-TS.Flux[Epnt-1][Tpnt-1]; 818 R-=TS.Flux[Epnt-1][Tpnt-1]; 819 820 R/=Tspread; 821 *TV+=TRange*R; 822 823 824 R1=TS.EInt[Epnt-1]+Espread*rand01(); 825 iStart=Epnt>lowPoly ? Epnt-lowPoly : 0; // max(Epnt-halfPoly,0) 826 iEnd=TS.nEnergy>Epnt+highPoly ? Epnt+highPoly : TS.nEnergy-1; // min(nEnergy-1,Epnt+highPoly 827 828 *EV=polInterp(TS.EInt+iStart,TS.EnergyBin+iStart,1+iEnd-iStart,R1); 829 830 if(*TV < TS.TimeBin[Tpnt-1] || *TV > TS.TimeBin[Tpnt]) 831 { 832 fprintf(stderr,"%d Tpnt %d Tval %g Epnt %d \n",TS.nTime,Tpnt,*TV,Epnt); 833 fprintf(stderr,"TBoundary == %12.6e,%g , %12.6e \n\n",TS.TimeBin[Tpnt-1],*TV,TS.TimeBin[Tpnt]); 834 } 835 836 837 if(*EV < *lim1 || *EV > *lim2) 838 { 839 fprintf(stderr,"outside boundaries\n Epnt= %d, Tpnt= %d binlo %g|%g| binhi %g \n", 840 Epnt,Tpnt,TS.EnergyBin[Epnt-1],*EV,TS.EnergyBin[Epnt]); 841 842 fprintf(stderr,"TS == %g %g :: %d %d \n",TS.EInt[Epnt-1],TS.EInt[Epnt],iStart,iEnd); 843 fprintf(stderr,"Points (%g) == ",R1); 844 for(i=0;i<iEnd-iStart;i++) 845 fprintf(stderr,"Time start %g %g",*(TS.EInt+i+iStart),*(TS.EnergyBin+iStart+i)); 846 fprintf(stderr,"\n"); 847 exit(1); 848 } 849 return; 850} 851 852int 853cmdnumberI(char *mc,int* num,const int len) 854 /*! 855 \param mc == character string to use 856 \param num :: Place to put output 857 \param len == length of the character string to process 858 returns 1 on success and 0 on failure 859 */ 860{ 861 int i,j; 862 char* ss; 863 char **endptr; 864 double nmb; 865 866 if (len<1) 867 return 0; 868 j=0; 869 870 for(i=0;i<len && mc[i] && 871 (mc[i]=='\t' || mc[i]==' ' || mc[i]==',');i++); 872 if(i==len || !mc[i]) return 0; 873 ss=malloc(sizeof(char)*(len+1)); 874 /* char *ss=new char[len+1]; */ 875 for(;i<len && mc[i]!='\n' && mc[i] 876 && mc[i]!='\t' && mc[i]!=' ' && mc[i]!=',';i++) 877 { 878 ss[j]=mc[i]; 879 j++; 880 } 881 if (!j) 882 { 883 free(ss); 884 return 0; //This should be impossible 885 } 886 ss[j]=0; 887 endptr=malloc(sizeof(char*)); 888 nmb = strtod(ss,endptr); 889 if (*endptr != ss+j) 890 { 891 free(endptr); 892 free(ss); 893 return 0; 894 } 895 *num = (double) nmb; 896 for(j=0;j<i && mc[j];j++) 897 mc[j]=' '; 898 free(endptr); 899 free(ss); 900 return 1; 901} 902 903 904FILE* 905openFile(char* FileName) 906/* 907 Tries to open the file: 908 (i) In current working directory 909 (ii) In MC_Path directory + ISIS_tables extension 910*/ 911{ 912 FILE* efile=0; 913 int fLen; 914 char extFileName[1024]; 915 char testFileName[2048]; 916 char mct[2048]; 917 918 fLen=strlen(FileName); 919 920 if (fLen>512) 921 { 922 fprintf(stderr,"Filename excessivley long [EXIT]:\n %s\n",FileName); 923 exit(1); 924 } 925 926 927 strcpy(extFileName,FileName); 928 strcpy(extFileName+fLen,".mcstas"); 929 930 /* Now check for the requested moderator file. In terms of precedence, 931 1) Use MCTABLES location if available 932 2) Is a local file available in PWD? 933 3) Is there an ISIS_tables in PWD? 934 4) Is the file available from the MCSTAS/contrib/ISIS_tables folder? 935 - otherwise fail */ 936 937 fprintf(stderr,"Searching for %s in multiple locations... -- \n",FileName); 938 939 /* 1) Is MCTABLES set and file located there? */ 940 if (getenv("MCTABLES")) 941 { 942 strcpy(mct, getenv("MCTABLES")); 943 sprintf(testFileName, "%s%c%s", mct, MC_PATHSEP_C, FileName); 944 efile=fopen(testFileName,"r"); 945 if (efile) { 946 fprintf(stderr," - Found in MCTABLES folder %s!\n",mct); 947 return efile; 948 } 949 } 950 951 /* 2) Is the file located in working dir? */ 952 efile=fopen(FileName,"r"); 953 if (efile) { 954 fprintf(stderr," - Found in current directory!\n"); 955 return efile; 956 } 957 958 efile=fopen(extFileName,"r"); 959 if (efile) return efile; 960 961 /* 3) Is the file in a local 'tablefiles' folder? */ 962 sprintf(testFileName,"%s%c%s","ISIS_tables",MC_PATHSEP_C,FileName); 963 efile=fopen(testFileName,"r"); 964 if (efile) { 965 fprintf(stderr," - Found in local ISIS_tables directory!\n"); 966 return efile; 967 } 968 969 /* 4) Is the file available within the MCSTAS install dir? */ 970 sprintf(testFileName,"%s%c%s%c%s%c%s",MCSTAS,MC_PATHSEP_C,"contrib",MC_PATHSEP_C,"ISIS_tables",MC_PATHSEP_C,FileName); 971 efile=fopen(testFileName,"r"); 972 if (efile) { 973 fprintf(stderr," - Found in MCSTAS system dir: \n %s%c%s%c%s%\n",MCSTAS,MC_PATHSEP_C,"contrib",MC_PATHSEP_C,"ISIS_tables"); 974 return efile; 975 } 976 977 /* If we reach here, the file was not found - raise fatal error */ 978 fprintf(stderr,"% - Not found! \nPlease check your McStas installation and/or MCTABLES environment variable!\n",FileName); 979 exit(1); 980} 981 982double strArea() 983{ 984 /* 985 Returns the mean Str view of the viewport 986 This integrates over each point on the window focus_xw to focus_yh 987 View port is symmetric so use only 1/4 of the view 988 for the calcuation. 989 Control Values TS.XAxis TS.ZAxis focus_xw focus_yh 990 */ 991 992 double A; 993 double Vx,Vy; // view temp points 994 double Mx,My; // moderator x,y 995 double D2; // Distance ^2 996 double projArea; // viewport projection to moderator 997 int i,j,aa,bb; // loop variables 998 int NStep; 999 1000 NStep=50; 1001 D2=dist*dist; 1002 A=0.0; 1003 fprintf(stderr,"TS axis == %g %g\n",TS.XAxis,TS.ZAxis); 1004 fprintf(stderr,"AW axis == %g %g\n",focus_xw,focus_yh); 1005 for(i=0;i<NStep;i++) // Mod X 1006 { 1007 Mx=(i+0.5)*TS.XAxis/(NStep*2.0); 1008 for(j=0;j<NStep;j++) // Mod Y 1009 { 1010 My=(j+0.5)*TS.ZAxis/(NStep*2.0); 1011 // Position on moderator == (Mx,My) 1012 for(aa= -NStep;aa<=NStep;aa++) //view port 1013 for(bb= -NStep;bb<=NStep;bb++) 1014 { 1015 Vx=aa*focus_xw/(2.0*NStep+1.0); 1016 Vy=bb*focus_yh/(2.0*NStep+1.0); 1017 A+=1.0/((Mx-Vx)*(Mx-Vx)+(My-Vy)*(My-Vy)+D2); 1018 } 1019 } 1020 } 1021 //change to Mx*My 1022 A*= (TS.XAxis*TS.ZAxis)/(NStep*NStep*(2.0*NStep+1.0)*(2.0*NStep+1.0)); 1023 // Correct for the area of the viewport. (tables are per cm^2) 1024 A*=focus_xw*focus_yh*10000; 1025 1026 // testing only - Goran 1027 // if (!modPosition) 1028 // { 1029 // projArea=focus_xw*focus_yh/tally/tally*(tally+dist)*(tally+dist); 1030 // A*=TS.XAxis*TS.ZAxis/(TS.XAxis*TS.ZAxis-projArea); 1031 // } 1032 1033 fprintf(stderr,"Viewport == %g %g Moderator size == (%g * %g) m^2 \n",focus_xw,focus_yh,TS.XAxis,TS.ZAxis); 1034 fprintf(stderr,"Dist == %g (metres) \n",dist); 1035 fprintf(stderr,"Viewport Solid angle == %g str\n",A/(focus_xw*focus_yh*10000.0)); 1036 fprintf(stderr,"Solid angle used == %g str\n",A); 1037 return A; 1038} 1039 1040%} 1041 1042INITIALIZE 1043%{ 1044 /* READ IN THE ENERGY FILE */ 1045 FILE* Tfile; 1046 1047 1048 Nsim=mcget_ncount(); // Number of points requested. 1049 1050 Tfile=openFile(Face); // Get open file 1051 rtE0=convertEnergy(E0); 1052 rtE1=convertEnergy(E1); 1053 orderEnergy(&rtE0,&rtE1); 1054 1055 readHtable(Tfile,rtE0,rtE1); 1056 fclose(Tfile); 1057 1058 /**********************************************************************/ 1059 1060 Tnpts=0; 1061 Ncount=0; 1062 1063 fprintf(stderr,"Face == %s \n",Face); 1064 fprintf(stderr,"Number of Energy Points == %d\n",TS.nEnergy); 1065 if (dist<0) 1066 { 1067 dist=TS.rdumMid; 1068 fprintf(stderr,"Setting distance to moderator surface == %g\n", 1069 dist); 1070 } 1071 /* Do solid angle correction */ 1072 angleArea= strArea(); 1073 1074 fprintf(stderr,"Totals:: %g %d %d \n",TS.Total,TS.nEnergy,TS.nTime); 1075 1076%} 1077 1078TRACE 1079%{ 1080 double v,r,E; 1081 double xf,yf,dx,dy; /* mxp ->max var in param space */ 1082 double Ival,Tval,Eval; 1083 1084 Ncount++; 1085 1086 x += TS.XAxis*(0.5-rand01()); /* Get point on shutter enterance */ 1087 y += TS.ZAxis*(0.5-rand01()); /* Get point on shutter enterance */ 1088 1089 // if (modPosition) z=TS.rdumMid; // Goran 1090 1091 xf = focus_xw*(0.5-rand01()); /* Choose focusing position uniformly */ 1092 yf = focus_yh*(0.5-rand01()); /* Choose focusing position uniformly */ 1093 1094 getPoint(&Tval,&Eval,&rtE0,&rtE1); 1095 1096// Ival=TS.Total*6.2415093e+12; // * 1.1879451; /* ( of proton in 1uAmp) * (1-cos(30))*2*Pi */ 1097// Ival=TS.Total*5.9294338e+12; // Number of proton in 1uAmp * 0.95 (0.95 because of muon target loss) 1098 Ival=TS.Total*6.2415093e+12; // Number of proton in 1uAmp 1099 1100 dx = xf-x; 1101 dy = yf-y; 1102 r = sqrt(dx*dx+dy*dy+dist*dist); // Actual distance to point 1103 v = SE2V*sqrt(Eval); // Calculate the velocity 1104 vz = v*fabs(dist)/r; 1105 vy = v*dy/r; 1106 vx = v*dx/r; 1107 1108 1109 t=Tval-(TS.rdumMid-TS.timeOffset)/vz; ///vz; 1110 1111 if (modPosition) 1112 { 1113 t+=TS.rdumMid/vz; 1114 } 1115 1116 p=beamcurrent*angleArea*Ival/Nsim; 1117 1118 // testing only - Goran 1119 // if (!modPosition) p*=TS.timeOffset; 1120 1121 1122 if (verbose && !(Ncount % 100000)) 1123 { 1124 fprintf(stderr,"FPos[%d]=> %g %g %g \n",Ncount,x,y,z); 1125 fprintf(stderr,"FDir[%d]=> %g %g %g \n",Ncount,vx,vy,vz); 1126 fprintf(stderr,"PlaneAxis %g %g \n",TS.XAxis,fullAngle); 1127 fprintf(stderr,"RMID %g \n",TS.rdumMid); 1128 fprintf(stderr,"TimeMid[%d]=> %g\n",Ncount,TS.rdumMid); 1129 fprintf(stderr,"TimeOffset[%d]=> %g\n",Ncount,TS.timeOffset); 1130 fprintf(stderr,"TimeTval[%d]=> %g\n",Ncount,Tval); 1131 fprintf(stderr,"TimeZero[%d]=> %g\n",Ncount,t); 1132 } 1133 1134%} 1135 1136MCDISPLAY 1137%{ 1138 double cirp=0.0,cirq=0.3,pi=3.141592654; 1139 int pp=0; /* circle drawing parameter*/ 1140 1141 1142 1143 magnify("xy"); 1144 multiline(5,-0.5*TS.XAxis,-0.5*TS.ZAxis,0.0, 1145 0.5*TS.XAxis,-0.5*TS.ZAxis,0.0, 1146 0.5*TS.XAxis,0.5*TS.ZAxis,0.0, 1147 -0.5*TS.XAxis,0.5*TS.ZAxis,0.0, 1148 -0.5*TS.XAxis,-0.5*TS.ZAxis,0.0); 1149 /* circle("xy",0.0,0.0,0.0,cos(cirp)); */ 1150 1151 /*line(0.5*sin(cirp),0.0,0.5*cos(cirp),0.5*sin(cirq),0.0,0.5*cos(cirq));*/ 1152 1153 /*line(-0.5,0.0,0.0,0.0,0.0,0.5);*/ 1154 1155 for (pp=0;pp<=20;pp=pp+2) 1156 { 1157 cirp= (pp*(pi/21.0))-(0.5*pi); 1158 cirq= ((pp+1)*(pi/21.0))-(0.5*pi); 1159 line(0.5*sin(cirp),0.0,0.5*cos(cirp),0.5*sin(cirq),0.0,0.5*cos(cirq)); 1160 } 1161 1162 %} 1163 1164END 1165