1/******************************************************************************* 2* 3* McStas, neutron ray-tracing package 4* Copyright 1997-2012, All rights reserved 5* Risoe National Laboratory, Roskilde, Denmark 6* Institut Laue Langevin, Grenoble, France 7* 8* Component: TOF_PSD_monitor_q 9* 10* %I 11* Based on: Henrich Frielinghaus, FZJuelich 12* Date: Apr 2013 13* Origin: xxx 14* Modified by: xxx 15* 16* Multiple TOF detectors for SANS instrument. 17* The component is to be placed at the sample position. 18* For the time being better switch gravity off. 19* 20* %D 21* TOF monitor that calculates I of q. 22* 23* Example: TOFSANSdet(); 24* 25* %P 26* INPUT PARAMETERS: 27* plenght: [s] pulse lenght 28* ssdist: [] Source-Sample distance for TOF calculation. 29* coldis: [] Collimation length 30* Sthckn: [cm] sample thickness 31* ds1: [] distance of detector 1 32* xw1: [] width of detector 1 (0 for off) 33* yh1: [] height of detector 1 (0 for off) 34* hl1: [] w/h of hole in det 1 (0 for off), full width 35* ds2: [] distance...........2 36* xw2: [] width..............2 37* yh2: [] heigth.............2 38* hl2: [] hole...............2 39* ds1: [] distance...........3 40* xw3: [] width..............3 41* yh3: [] height.............3 42* hl3: [] hole...............3 (beam stop, used for primary beam detection) 43* vx3: [] vertical extension of beam stop down (in case of gravity off set 0.0, otherwise value larger than 1.0) 44* tmin: [s] Beginning of time window 45* tmax: [s] End of time window 46* Nx: [] Number of horizontal detector pixels (detector 1,2,3) better leave unchanged 47* Ny: [] Number of vertical detector pixels (detector 1,2,3) better leave unchanged 48* Nt: [] Number of time bins (detector 1,2,3) better leave unchanged 49* qmin: [1/A] Lower limit of q-range 50* qmax: [1/A] Upper limit of q-range 51* Nq: [1] Number of q-bins 52* fname: [] file name (first part without extensions) 53* rstneu: restore neutron after treatment ??? (0.0 = no) 54* centol: [] tolerance of center determination (if center larger than centol*calculated_center then set back to theory) 55* inttol: [] tolerance of intensity (if primary beam intensity smaller than inttol*max_intensity then discard data) 56* qcal: [] calibration of intensity (cps) to width of q-bin (non_zero = yes, zero = no) 57* 58* OUTPUT PARAMETERS: 59* 60* PSDq_N: [] Array of neutron counts 61* PSDq_p: [] Array of neutron weight counts 62* PSDq_p2: [] Array of second moments 63* 64* %E 65*******************************************************************************/ 66 67DEFINE COMPONENT TOFSANSdet 68DEFINITION PARAMETERS (Nq=100, string fname) 69SETTING PARAMETERS (plength=0.00286, ssdist=27.0, coldis = 20.0, Sthckn = 0.1, 70ds1=1.0, xw1=0.0, yh1=0.0, hl1=0.0, /* first one switched off */ 71ds2=5.0, xw2=1.0, yh2=1.0, hl2=0.2, 72ds3=20.0, xw3=1.0, yh3=1.0, hl3=0.04, 73vx3=0.0, 74tmin=0.005, tmax=0.15, 75Nx=128.0, Ny=128.0, Nt=500.0, /* be cautious with large numbers */ 76qmin=0.0005, qmax=0.11, // Nq=500.0, <<<<<<<<<<<<<<<<< defined above 77rstneu = 0.0, 78centol = 0.1, inttol = 0.0001, qcal=1.0) 79OUTPUT PARAMETERS (simplN,simplI,simplIe) 80 81SHARE 82%{ 83#include <stdlib.h> 84#define TSdN(i,j,s,k) TSdNf[i+Nxc*(j+Nyc*(s+Ntc*(k)))] 85#define TSdp(i,j,s,k) TSdpf[i+Nxc*(j+Nyc*(s+Ntc*(k)))] 86#define TSdp2(i,j,s,k) TSdp2f[i+Nxc*(j+Nyc*(s+Ntc*(k)))] 87#define calibQ(i,j) calibQf[i+Nqc*(j)] 88#define calibN(i,j) calibNf[i+Nqc*(j)] 89#define calibI(i,j) calibIf[i+Nqc*(j)] 90#define calibIe(i,j) calibIef[i+Nqc*(j)] 91#define calibWt(i,j) calibWtf[i+Nqc*(j)] 92#define sectnQ(i,j) sectnQf[i+Nqc*(j)] 93#define sectnN(i,j) sectnNf[i+Nqc*(j)] 94#define sectnI(i,j) sectnIf[i+Nqc*(j)] 95#define sectnIe(i,j) sectnIef[i+Nqc*(j)] 96#define sectnWt(i,j) sectnWtf[i+Nqc*(j)] 97 98double ddmin(double A, double B) { 99if (A<B) return A; else return B; 100}; 101 102double ddmax(double A, double B) { 103if (A>B) return A; else return B; 104}; 105 106int iimin(int A, int B) { 107if (A<B) return A; else return B; 108}; 109 110int iimax(int A, int B) { 111if (A>B) return A; else return B; 112}; 113 114double storez,storen; 115 116%} 117 118DECLARE 119%{ 120 int Nxc,Nyc,Ntc,Nqc; 121 double *TSdNf; 122 double *TSdpf; 123 double *TSdp2f; 124 int i,j,s,k; 125 double ds1c,ds2c,ds3c; 126 double xw1c,xw2c,xw3c; 127 double yh1c,yh2c,yh3c; 128 double hl1c,hl2c,hl3c; 129 double vx3c; 130 double Pic; 131 132// double simplQ[Nq]; 133 double simplN[Nq]; 134 double simplI[Nq]; 135 double simplIe[Nq]; 136 137 double Ncount,Pcount; 138%} 139 140 141INITIALIZE 142%{ 143 storez = 0.0; 144 storen = 0.0; 145 146 Nxc = floor( ddmin(ddmax(Nx,16.0),512.0) + 0.50); 147 Nyc = floor( ddmin(ddmax(Ny,16.0),512.0) + 0.50); 148 Ntc = floor( ddmin(ddmax(Nt,10.0),2000.0)+ 0.50); 149 Nqc = floor( ddmin(ddmax(Nq,1.00),2000.0)+ 0.50); 150 151 TSdNf = (double*)malloc(sizeof(double)*Nxc*Nyc*Ntc*3); 152 TSdpf = (double*)malloc(sizeof(double)*Nxc*Nyc*Ntc*3); 153 TSdp2f= (double*)malloc(sizeof(double)*Nxc*Nyc*Ntc*3); 154 155 for (i=0; i<Nxc; i++) 156 for (j=0; j<Nyc; j++) 157 for (s=0; s<Ntc; s++) 158 for(k=0; k<3; k++) 159 { 160 TSdN(i,j,s,k) = 0.0; 161 TSdp(i,j,s,k) = 0.0; 162 TSdp2(i,j,s,k) = 0.0; 163 }; 164 165 ds1c = ddmax(ds1,0.0); 166 ds2c = ddmax(ds2,0.0); 167 ds3c = ddmax(ds3,0.0); 168 if (ds2c>=ds3c) ds2c = 0.0; 169 if (ds1c>=ds2c) ds1c = 0.0; 170 xw1c = ddmax(xw1,0.0); 171 yh1c = ddmax(yh1,0.0); 172 xw2c = ddmax(xw2,0.0); 173 yh2c = ddmax(yh2,0.0); 174 xw3c = ddmax(xw3,0.0); 175 yh3c = ddmax(yh3,0.0); 176 hl1c = hl1; 177 hl2c = hl2; 178 hl3c = hl3; 179 if (ds3c==0.0) {xw3c=0.0; yh3c=0.0; hl3c=0.0;}; 180 if (ds2c==0.0) {xw2c=0.0; yh2c=0.0; hl2c=0.0;}; 181 if (ds1c==0.0) {xw1c=0.0; yh1c=0.0; hl1c=0.0;}; 182 if (hl1c<=0.0) {hl1c = 0.0;} else {hl1c=ddmax(hl1c,ddmax(xw1c/Nx,yh1c/Ny)*3.0);}; 183 if (hl2c<=0.0) {hl2c = 0.0;} else {hl2c=ddmax(hl2c,ddmax(xw2c/Nx,yh2c/Ny)*3.0);}; 184 hl3c=ddmax(hl3c,ddmax(xw3c/Nx,yh3c/Ny)*3.0); /* leave some space for other data */ 185 vx3c = ddmin(ddmax(vx3,1.0),15.0); 186 187/* if (fname.empty()) fname="SANSareaDet"; */ 188 189 Pic = 3.141592653589793238462643; 190 191 Ncount = 0.0; 192 Pcount = 0.0; 193%} 194 195TRACE 196 %{ 197 int i,j,s,k; 198 double tt, zpos; 199 double absflg; 200 201 PROP_Z0; 202 zpos = 0.0; 203 absflg = 0.0; 204 205 storez += vz*t; 206 storen++; 207 208// if (p<0.0) fprintf(stdout,"x"); 209 210 k = 0; 211 if (xw1c>0.0 && yh1c>0.0) { 212 tt = (ds1c-zpos)/vz; // for all functionality the gravity direction should be in -y direction 213 PROP_DT(tt); 214 zpos = ds1c; 215 if (fabs(x)<0.5*xw1c && fabs(y)<0.5*yh1c && (fabs(x)>0.5*hl1c || fabs(y)>0.5*hl1c)) { 216 tt = t - 0.5*plength; // Actual time of flight minus one half pulsewidth. 217 absflg = 1.0; 218 s = floor( (tt-tmin) *Ntc/(tmax-tmin) ); /* Bin number */ 219 i = floor( (x+0.5*xw1c)*Nxc/ xw1c ); 220 j = floor( (y+0.5*yh1c)*Nyc/ yh1c ); 221 if (s>=0 && s<Ntc) { 222 TSdN(i,j,s,k)++; 223 TSdp(i,j,s,k) += p; 224 TSdp2(i,j,s,k)+= p*p; 225 }; 226 SCATTER; 227 }; 228 }; 229 230 k = 1; 231 if (xw2c>0.0 && yh2c>0.0 && absflg==0.0) { 232 tt = (ds2c-zpos)/vz; 233 PROP_DT(tt); 234 zpos = ds2c; 235 if (fabs(x)<0.5*xw2c && fabs(y)<0.5*yh2c && (fabs(x)>0.5*hl2c || fabs(y)>0.5*hl2c)) { 236 tt = t - 0.5*plength; // Actual time of flight minus one half pulsewidth. 237 absflg = 1.0; 238 s = floor( (tt-tmin) *Ntc/(tmax-tmin) ); /* Bin number */ 239 i = floor( (x+0.5*xw2c)*Nxc/ xw2c ); 240 j = floor( (y+0.5*yh2c)*Nyc/ yh2c ); 241 if (s>=0 && s<Ntc) { 242 TSdN(i,j,s,k)++; 243 TSdp(i,j,s,k) += p; 244 TSdp2(i,j,s,k)+= p*p; 245 }; 246 SCATTER; 247 }; 248 }; 249 250 k = 2; 251 if (xw3c>0.0 && yh3c>0.0 && absflg==0.0) { 252 tt = (ds3c-zpos)/vz; 253 PROP_DT(tt); 254 zpos = ds3c; 255 if (fabs(x)<0.5*xw3c && fabs(y)<0.5*yh3c && (fabs(x)>0.5*hl3c || y>0.5*hl3c || y<(0.5-vx3c)*hl3c)) { 256 tt = t - 0.5*plength; // Actual time of flight minus one half pulsewidth. 257 absflg = 1.0; 258 s = floor( (tt-tmin) *Ntc/(tmax-tmin) ); /* Bin number */ 259 i = floor( (x+0.5*xw3c)*Nxc/ xw3c ); 260 j = floor( (y+0.5*yh3c)*Nyc/ yh3c ); 261 if (s>=0 && s<Ntc) { 262 TSdN(i,j,s,k)++; 263 TSdp(i,j,s,k) += p; 264 TSdp2(i,j,s,k)+= p*p; 265 }; 266 SCATTER; 267 }; 268 }; 269 270 Ncount += absflg; 271 272 k = 2; 273 tt = (ds3c-zpos)/vz; 274 if (tt>0.0) PROP_DT(tt); 275 zpos = ds3c; 276 if (fabs(x)<=0.5*hl3c && y<=0.5*hl3c && y>=(0.5-vx3c)*hl3c) { 277 tt = t - 0.5*plength; // Actual time of flight minus one half pulsewidth. 278 absflg = 1.0; 279 s = floor( (tt-tmin) *Ntc/(tmax-tmin) ); /* Bin number */ 280 i = Nxc/2; 281 j = Nyc/2; 282 if (s>=0 && s<Ntc) { 283 TSdN(i,j,s,k)++; 284 TSdp(i,j,s,k) += p; 285 TSdp2(i,j,s,k)+= p*p; 286 TSdN(i-1,j,s,k)++; /* weight for x */ 287 TSdp(i-1,j,s,k) += p; 288 TSdp2(i-1,j,s,k)+= p*x; 289 TSdN(i,j-1,s,k)++; /* weight for y */ 290 TSdp(i,j-1,s,k) += p; 291 TSdp2(i,j-1,s,k)+= p*y; 292 }; 293 SCATTER; 294 Pcount++; 295 }; 296 297/* if (absflg!=0.0) ABSORB; */ 298 299 if (rstneu!=0.0) { RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); }; 300 301 %} 302 303SAVE 304 %{ 305 int i,j,s,k,qi,si; 306// int kkk; 307 double dsA[3],xwA[3],yhA[3],hlA[3]; 308 dsA[0] = ds1c; 309 dsA[1] = ds2c; 310 dsA[2] = ds3c; 311 xwA[0] = xw1c; 312 xwA[1] = xw2c; 313 xwA[2] = xw3c; 314 yhA[0] = yh1c; 315 yhA[1] = yh2c; 316 yhA[2] = yh3c; 317 hlA[0] = hl1c; 318 hlA[1] = hl2c; 319 hlA[2] = hl3c; 320 321 double q1,q2,qstp,qminl,qmaxl,qcalf; 322 double time,time3,s3cl,s3sw; 323 int s3,s3sg,kfl; 324 double wght1,wght2; 325 double xcen,ycen,prIn,vvz,yguess,maxint; 326 double ic1,ic2,ic3,ic4,jc1,jc2,jc3,jc4; 327 int imn,imx,im3,im4,jmn,jmx,jm3,jm4; 328 double cimn,cimx,cjmn,cjmx; 329 double Xi,Yj,Ni,Nj,delN,scali,scalj,scal; 330 double TNw2, CSw, SNwh; 331 double dsA2, DIS1, DIS2, LAM, QQQ, OMG, AREA, AR2, FAK, determ; 332 double Q11,Q21,Q12,Q22; 333 double N11,N21,N12,N22; 334 double Qx,Qy; 335 double Nmn,Nmx,T00,Tmn,Tmx; 336 double Nit,Njt,NNt,Ttt,Qstp; 337 double Qmn,Qmx; 338 int qq1,qq2; 339 double qmn,qmx,qst,qqq; 340 double tmn,tmx,tst,ttt; 341 double NNN,Int,Err,Wt; 342 double Qxx,Qyy,Qdx,Qdy; 343 344 char filename[99]; 345 346// double* calibQf; 347 double* calibNf; 348 double* calibIf; 349 double* calibIef; 350 double* calibWtf; 351// double* sectnQf; 352 double* sectnNf; 353 double* sectnIf; 354 double* sectnIef; 355 double* sectnWtf; 356 357// simplQ = (double*)malloc(sizeof(double)*Nqc); 358// simplN = (double*)malloc(sizeof(double)*Nqc); 359// simplI = (double*)malloc(sizeof(double)*Nqc); 360// simplIe = (double*)malloc(sizeof(double)*Nqc); 361// calibQf = (double*)malloc(sizeof(double)*Nqc*4); 362 calibNf = (double*)malloc(sizeof(double)*Nqc*4); 363 calibIf = (double*)malloc(sizeof(double)*Nqc*4); 364 calibIef= (double*)malloc(sizeof(double)*Nqc*4); 365 calibWtf= (double*)malloc(sizeof(double)*Nqc*4); 366// sectnQf = (double*)malloc(sizeof(double)*Nqc*20); 367 sectnNf = (double*)malloc(sizeof(double)*Nqc*20); 368 sectnIf = (double*)malloc(sizeof(double)*Nqc*20); 369 sectnIef= (double*)malloc(sizeof(double)*Nqc*20); 370 sectnWtf= (double*)malloc(sizeof(double)*Nqc*20); 371 372 storez /= storen; 373// if (ssdist>0.1*storez) storez = ssdist; 374 375// fprintf(stdout,"%g %g \n",Ncount,Pcount); 376// fprintf(stdout,"point1"); 377// fprintf(stdout,"xxxx %g %g \n",storez,ssdist); 378 379 i = Nxc/2; /* calculate center values on beam stop */ 380 j = Nyc/2; 381 k = 2; 382 maxint = 0.0; 383 for (s=0; s<Ntc; s++) { 384 if (TSdp(i,j,s,k)>0.0) { 385 TSdp2(i-1,j,s,k) /= TSdp(i-1,j,s,k); 386 TSdp2(i,j-1,s,k) /= TSdp(i,j-1,s,k); 387 maxint = ddmax(maxint,TSdp(i,j,s,k)); 388// fprintf(stdout," %g ",maxint); 389 }; 390 }; 391 maxint *= inttol; 392// fprintf(stdout," %g ",maxint); 393 394// fprintf(stdout,"point2"); 395 396 qstp = log(qmax/qmin)/Nqc; 397 qminl= log(qmin); 398 qmaxl= log(qmax); 399 400 for (qi=0; qi<Nqc; qi++) { 401// q1 = qmin * exp(qi*qstp); 402// q2 = q1 * exp(qstp); 403// simplQ[qi] = 0.5*(q1+q2); 404 simplN[qi] = 0.0; 405 simplI[qi] = 0.0; 406 simplIe[qi]= 0.0; 407 for (si=0; si<4; si++) { 408// calibQ(qi,si) = simplQ[qi]; 409 calibN(qi,si) = 0.0; 410 calibI(qi,si) = 0.0; 411 calibIe(qi,si)= 0.0; 412 calibWt(qi,si)= 0.0; 413 }; 414 for (si=0; si<20; si++) { 415// sectnQ(qi,si) = simplQ[qi]; 416 sectnN(qi,si) = 0.0; 417 sectnI(qi,si) = 0.0; 418 sectnIe(qi,si)= 0.0; 419 sectnWt(qi,si)= 0.0; 420 }; 421 }; 422 423// fprintf(stdout,"point3"); 424 425 kfl = 0; 426 for (k=0; k<3; k++) { 427// fprintf(stdout,"pointK"); 428 if (dsA[k]>0.0 && xwA[k]>0.0 && yhA[k]>0.0) { 429 for (s=0; s<Ntc; s++) { 430// fprintf(stdout,"S"); 431 time = tmin + (s+0.5)*(tmax-tmin)/Ntc; 432 time3 = time*(ssdist+ds3c)/(ssdist+dsA[k]); 433 s3cl = (time3-tmin)*Ntc/(tmax-tmin); 434 s3 = floor(s3cl); 435 s3sw = s3cl-s3-0.5; 436 if (s3sw<0.0) {s3sg = -1;} else {s3sg = 1;}; 437 if (s3>=0 && s3<Ntc) { 438 i = Nxc/2; 439 j = Nyc/2; 440 xcen = TSdp2(i-1,j,s3,2); 441 ycen = TSdp2(i,j-1,s3,2); 442 prIn = TSdp(i,j,s3,2); 443 if (prIn>=maxint) {wght1=1.0-fabs(s3sw);} else {wght1=0.0;}; 444 wght2 = 0.0; 445 if (s3+s3sg>=0 && s3+s3sg<Ntc) { 446 if (TSdp(i,j,s3+s3sg,2)>=maxint) wght2=fabs(s3sw); 447 if (wght1+wght2>0.0) { 448 xcen = (wght2*TSdp2(i-1,j,s3+s3sg,2) + wght1*xcen)/(wght1+wght2); 449 ycen = (wght2*TSdp2(i,j-1,s3+s3sg,2) + wght1*ycen)/(wght1+wght2); 450 prIn = (wght2*TSdp(i,j,s3+s3sg,2) + wght1*prIn)/(wght1+wght2); 451 } else { 452 xcen = 0.0; 453 ycen = 0.0; 454 prIn = 0.0; 455 }; 456 }; 457 prIn *= (ssdist+ds3c)/(ssdist+dsA[k]); // correct for spreading of times 458// fprintf(stdout,"%g %g \n",maxint,prIn); 459 if (fabs(xcen)>hlA[2]*centol) {xcen=0.0;} else {xcen *= dsA[k]/ds3c;}; 460 if (mcgravitation==0) {if (fabs(ycen)>hlA[2]*centol) {ycen=0.0;} else {ycen *= dsA[k]/ds3c;}; } 461 else {vvz = (ssdist+ds3c)/time3; 462 yguess = 0.5*GRAVITY*(pow(0.5*coldis+ds3c,2)-0.25*coldis*coldis)/(vvz*vvz); 463 if (fabs(ycen-yguess)>hlA[2]*centol) {ycen=yguess;} else 464 {ycen*=(pow(0.5*coldis+dsA[k],2)-0.25*coldis*coldis)/(pow(0.5*coldis+ds3c,2)-0.25*coldis*coldis);}; }; 465 if (kfl==0) { 466 ic1 = -0.5; 467 ic2 = Nxc-0.5; 468 jc1 = -0.5; 469 jc2 = Nyc-0.5; 470 imn = 0; 471 imx = Nxc-1; 472 jmn = 0; 473 jmx = Nyc-1; 474 } else { 475 ic1 = ddmax((-0.5*hlA[k-1]*dsA[k]/dsA[k-1]+xcen)*Nxc/xwA[k]+0.5*(Nxc-1), -0.5); 476 ic2 = ddmin(( 0.5*hlA[k-1]*dsA[k]/dsA[k-1]+xcen)*Nxc/xwA[k]+0.5*(Nxc-1),Nxc-0.5); 477 jc1 = ddmax((-0.5*hlA[k-1]*dsA[k]/dsA[k-1]+ycen)*Nyc/yhA[k]+0.5*(Nyc-1), -0.5); 478 jc2 = ddmin(( 0.5*hlA[k-1]*dsA[k]/dsA[k-1]+ycen)*Nyc/yhA[k]+0.5*(Nyc-1),Nyc-0.5); 479 imn = floor( ic1+0.5); 480 imx = -floor(-ic2+0.5); 481 jmn = floor( jc1+0.5); 482 jmx = -floor(-jc2+0.5); 483 }; 484 ic3 = -0.5*hlA[k]*Nxc/xwA[k]+0.5*(Nxc-1); 485 ic4 = 0.5*hlA[k]*Nxc/xwA[k]+0.5*(Nxc-1); 486 jc3 = -0.5*hlA[k]*Nyc/yhA[k]+0.5*(Nyc-1); 487 jc4 = 0.5*hlA[k]*Nyc/yhA[k]+0.5*(Nyc-1); 488 im3 = -floor(-ic3-0.5); 489 im4 = floor( ic4-0.5); 490 jm3 = -floor(-jc3-0.5); 491 jm4 = floor( jc4-0.5); 492 493// fprintf(stdout,"s3 %g %g %g %g %3i %3i %3i %3i \n",ic1,ic2,ic3,ic4,imn,imx,im3,im4); 494 495// if (imn<0 || jmn<0 || imx>Nxc-1 || jmx>Nyc-1) fprintf(stdout,"error in index imn..."); 496 497 for (i=imn; i<=imx; i++) 498 for (j=jmn; j<=jmx; j++) 499 { if (i<im3 || i>im4 || j<jm3 || j>jm4) { 500 501// if (i<5 && j<5) fprintf(stdout,"i"); 502 503 cimn = i-0.5; 504 cimx = i+0.5; 505 cjmn = j-0.5; 506 cjmx = j+0.5; 507 if (i==imn) cimn=ic1; 508 if (i==imx) cimx=ic2; 509 if (cimn<=ic3 && cimx>=ic3) cimx = ic3; 510 if (cimn<=ic4 && cimx>=ic4) cimn = ic4; 511 if (j==jmn) cjmn=jc1; 512 if (j==jmx) cjmx=jc2; 513 if (cjmn<=jc3 && cjmx>=jc3) cjmx = jc3; 514 if (cjmn<=jc4 && cjmx>=jc4) cjmn = jc4; 515 AREA = (cimx-cimn)*(cjmx-cjmn); 516 517 dsA2 = dsA[k]*dsA[k]; 518 519 delN = 0.1; /* variation size for derivatives */ 520 521 Xi= (-0.5+(i+0.5+delN)/Nxc)*xwA[k]-xcen; /* derivative in X */ 522 Yj= (-0.5+(j+0.5 )/Nyc)*yhA[k]-ycen; 523 DIS1 = Xi*Xi+Yj*Yj; 524 DIS2 = dsA2+DIS1; 525 LAM = (2.0*Pic/V2K)*time/(ssdist+sqrt(DIS2)); 526 TNw2 = DIS1/dsA2; /* tan (theta) squared */ 527 CSw = 1.0/sqrt(1.0+TNw2); /* cos (theta) */ 528 SNwh = sqrt(0.5*(1.0-CSw)); /* sin (theta/2) */ 529 QQQ = 4.0*Pic*SNwh/LAM; 530 Q21 = QQQ/sqrt(DIS1); 531 Q11 = Xi*Q21; 532 Q21 *= Yj; 533 534 Xi= (-0.5+(i+0.5 )/Nxc)*xwA[k]-xcen; /* derivative in Y */ 535 Yj= (-0.5+(j+0.5+delN)/Nyc)*yhA[k]-ycen; 536 DIS1 = Xi*Xi+Yj*Yj; 537 DIS2 = dsA2+DIS1; 538 LAM = (2.0*Pic/V2K)*time/(ssdist+sqrt(DIS2)); 539 TNw2 = DIS1/dsA2; /* tan (theta) squared */ 540 CSw = 1.0/sqrt(1.0+TNw2); /* cos (theta) */ 541 SNwh = sqrt(0.5*(1.0-CSw)); /* sin (theta/2) */ 542 QQQ = 4.0*Pic*SNwh/LAM; 543 Q22 = QQQ/sqrt(DIS1); 544 Q12 = Xi*Q22; 545 Q22 *= Yj; 546 547 Xi= (-0.5+(i+0.5)/Nxc)*xwA[k]-xcen; /* main Q */ 548 Yj= (-0.5+(j+0.5)/Nyc)*yhA[k]-ycen; 549 DIS1 = Xi*Xi+Yj*Yj; 550 DIS2 = dsA2+DIS1; 551 LAM = (2.0*Pic/V2K)*time/(ssdist+sqrt(DIS2)); 552 TNw2 = DIS1/dsA2; /* tan (theta) squared */ 553 CSw = 1.0/sqrt(1.0+TNw2); /* cos (theta) */ 554 SNwh = sqrt(0.5*(1.0-CSw)); /* sin (theta/2) */ 555 QQQ = 4.0*Pic*SNwh/LAM; 556 Qy = QQQ/sqrt(DIS1); 557 Qx = Xi*Qy; 558 Qy *= Yj; 559 560 Q11 = (Q11-Qx)/delN; 561 Q21 = (Q21-Qy)/delN; 562 Q12 = (Q12-Qx)/delN; 563 Q22 = (Q22-Qy)/delN; 564 565 determ = Q11*Q22-Q21*Q12; 566 N11 = Q22/determ; 567 N21 = -Q21/determ; 568 N12 = -Q12/determ; 569 N22 = Q11/determ; 570 571 OMG = CSw*xwA[k]*yhA[k]/((Nxc*Nyc)*DIS2); 572 573 qi = floor(log(QQQ/qmin)/qstp); /* dump the original intensities */ 574 if (qi>=0 && qi<Nqc) { 575 qcalf = 1.0; 576 if (qcal!=0.0) { 577 q1 = qmin * exp(qi*qstp); 578 q2 = q1 * exp(qstp); 579 qcalf = q2-q1; 580 }; 581 simplN[qi] += TSdN(i,j,s,k); 582 simplI[qi] += TSdp(i,j,s,k) / qcalf; 583 simplIe[qi]+= TSdp2(i,j,s,k)/(qcalf*qcalf); 584 }; 585 586 if (prIn>=maxint) { 587 588 if (yhA[k]/Nyc>xwA[k]/Nxc) { 589 scali = 1.0; 590 scalj = yhA[k]*Nxc/(xwA[k]*Nyc); 591 scal = xwA[k]/Nxc; 592 } else { 593 scali = xwA[k]*Nyc/(yhA[k]*Nxc); 594 scalj = 1.0; 595 scal = yhA[k]/Nyc; 596 }; 597 598 Ni = i+0.5-(0.5+xcen/xwA[k])*Nxc*scali; /* (0,0) center */ 599 Nj = j+0.5-(0.5+ycen/yhA[k])*Nyc*scalj; 600 Nmn = Ni*Ni+Nj*Nj; /* find minimal and maximal Qs and phis */ 601 Nmx = Nmn; 602 T00 = atan2(Nj,Ni); 603 Tmn = T00; 604 Tmx = T00; 605 606 Nit = Ni+scali; /* (1,0) */ 607 Njt = Nj; 608 NNt = Nit*Nit+Njt*Njt; 609 Ttt = atan2(Njt,Nit); 610 if (fabs(Ttt+2.0*Pic-T00)<fabs(Ttt-T00)) Ttt+=2.0*Pic; 611 if (fabs(Ttt-2.0*Pic-T00)<fabs(Ttt-T00)) Ttt-=2.0*Pic; 612 Nmn = ddmin(Nmn,NNt); 613 Nmx = ddmax(Nmx,NNt); 614 Tmn = ddmin(Tmn,Ttt); 615 Tmx = ddmax(Tmx,Ttt); 616 617// Nit = Ni+scali; /* (1,1) */ 618 Njt = Nj+scalj; 619 NNt = Nit*Nit+Njt*Njt; 620 Ttt = atan2(Njt,Nit); 621 if (fabs(Ttt+2.0*Pic-T00)<fabs(Ttt-T00)) Ttt+=2.0*Pic; 622 if (fabs(Ttt-2.0*Pic-T00)<fabs(Ttt-T00)) Ttt-=2.0*Pic; 623 Nmn = ddmin(Nmn,NNt); 624 Nmx = ddmax(Nmx,NNt); 625 Tmn = ddmin(Tmn,Ttt); 626 Tmx = ddmax(Tmx,Ttt); 627 628 Nit = Ni; /* (0,1) */ 629// Njt = Nj+scalj; 630 NNt = Nit*Nit+Njt*Njt; 631 Ttt = atan2(Njt,Nit); 632 if (fabs(Ttt+2.0*Pic-T00)<fabs(Ttt-T00)) Ttt+=2.0*Pic; 633 if (fabs(Ttt-2.0*Pic-T00)<fabs(Ttt-T00)) Ttt-=2.0*Pic; 634 Nmn = ddmin(Nmn,NNt); 635 Nmx = ddmax(Nmx,NNt); 636 Tmn = ddmin(Tmn,Ttt); 637 Tmx = ddmax(Tmx,Ttt); 638 639 Nit = Ni-scali; /* (-1,1) */ 640// Njt = Nj+scalj; 641 NNt = Nit*Nit+Njt*Njt; 642 Ttt = atan2(Njt,Nit); 643 if (fabs(Ttt+2.0*Pic-T00)<fabs(Ttt-T00)) Ttt+=2.0*Pic; 644 if (fabs(Ttt-2.0*Pic-T00)<fabs(Ttt-T00)) Ttt-=2.0*Pic; 645 Nmn = ddmin(Nmn,NNt); 646 Nmx = ddmax(Nmx,NNt); 647 Tmn = ddmin(Tmn,Ttt); 648 Tmx = ddmax(Tmx,Ttt); 649 650// Nit = Ni-scali; /* (-1,0) */ 651 Njt = Nj; 652 NNt = Nit*Nit+Njt*Njt; 653 Ttt = atan2(Njt,Nit); 654 if (fabs(Ttt+2.0*Pic-T00)<fabs(Ttt-T00)) Ttt+=2.0*Pic; 655 if (fabs(Ttt-2.0*Pic-T00)<fabs(Ttt-T00)) Ttt-=2.0*Pic; 656 Nmn = ddmin(Nmn,NNt); 657 Nmx = ddmax(Nmx,NNt); 658 Tmn = ddmin(Tmn,Ttt); 659 Tmx = ddmax(Tmx,Ttt); 660 661// Nit = Ni-scali; /* (-1,-1) */ 662 Njt = Nj-scalj; 663 NNt = Nit*Nit+Njt*Njt; 664 Ttt = atan2(Njt,Nit); 665 if (fabs(Ttt+2.0*Pic-T00)<fabs(Ttt-T00)) Ttt+=2.0*Pic; 666 if (fabs(Ttt-2.0*Pic-T00)<fabs(Ttt-T00)) Ttt-=2.0*Pic; 667 Nmn = ddmin(Nmn,NNt); 668 Nmx = ddmax(Nmx,NNt); 669 Tmn = ddmin(Tmn,Ttt); 670 Tmx = ddmax(Tmx,Ttt); 671 672 Nit = Ni; /* (0,-1) */ 673// Njt = Nj-scalj; 674 NNt = Nit*Nit+Njt*Njt; 675 Ttt = atan2(Njt,Nit); 676 if (fabs(Ttt+2.0*Pic-T00)<fabs(Ttt-T00)) Ttt+=2.0*Pic; 677 if (fabs(Ttt-2.0*Pic-T00)<fabs(Ttt-T00)) Ttt-=2.0*Pic; 678 Nmn = ddmin(Nmn,NNt); 679 Nmx = ddmax(Nmx,NNt); 680 Tmn = ddmin(Tmn,Ttt); 681 Tmx = ddmax(Tmx,Ttt); 682 683 Nit = Ni+scali; /* (1,-1) */ 684// Njt = Nj-scalj; 685 NNt = Nit*Nit+Njt*Njt; 686 Ttt = atan2(Njt,Nit); 687 if (fabs(Ttt+2.0*Pic-T00)<fabs(Ttt-T00)) Ttt+=2.0*Pic; 688 if (fabs(Ttt-2.0*Pic-T00)<fabs(Ttt-T00)) Ttt-=2.0*Pic; 689 Nmn = ddmin(Nmn,NNt); 690 Nmx = ddmax(Nmx,NNt); 691 Tmn = ddmin(Tmn,Ttt); 692 Tmx = ddmax(Tmx,Ttt); 693 694 DIS1 = Nmn*scal*scal; /* minimum Q provides largest distance */ 695 DIS2 = dsA2+DIS1; 696 LAM = (2.0*Pic/V2K)*time/(ssdist+sqrt(DIS2)); 697 TNw2 = DIS1/dsA2; /* tan (theta) squared */ 698 CSw = 1.0/sqrt(1.0+TNw2); /* cos (theta) */ 699 SNwh = sqrt(0.5*(1.0-CSw)); /* sin (theta/2) */ 700 Qmn = 4.0*Pic*SNwh/LAM; 701 Qmx = QQQ*QQQ/Qmn; 702 703// Nmn = sqrt(Nmn); // not needed anymore 704// Nmx = sqrt(Nmx); 705 706 LAM = (2.0*Pic/V2K)*time/(ssdist+sqrt(dsA2)); 707 Qstp = 2.0*Pic*min(xwA[k]/Nxc,yhA[k]/Nyc)/(sqrt(dsA2)*LAM); 708 709 qq1 = iimax( floor(log(Qmn/qmin)/qstp-1.0), 0); 710 qq2 = iimin(-floor(log(qmin/Qmx)/qstp-1.0),Nqc-1); 711 712// if (j==30 && i==30) fprintf(stdout,"%g %g %g %g %i %i \n",QQQ,Qmn,Qmx,Qstp,qq1,qq2); 713// if (j==30 && i==30) fprintf(stdout,"%g %g %g %g %g \n",QQQ,Q11,Q12,Q21,Q22); 714// if (j==30 && i==30 && LAM>4.08 && LAM<4.30) fprintf(stdout,"%g %g %g %i %i %g %g %g %g %g %g \n",LAM,OMG,QQQ,qq1,qq2,scali,scalj,N11,N22,Qstp,prIn); 715// if (j==30 && i==30 && LAM>4.08 && LAM<4.30) fprintf(stdout," %g %g %g \n",Qmn,Qmx,AREA); 716 717// kkk =0; 718 719 for (qi=qq1; qi<=qq2; qi++) { 720// if (qi==qq1) fprintf(stdout,"q"); 721 q1 = qmin * exp(qi*qstp); 722 q2 = q1 * exp(qstp); 723 if (q2-q1<=1.5*Qstp) { 724 qmn = sqrt(q1*q2); 725 qmx = qmn; 726 qst = Qstp; 727 qmx+= 1e-6*qst; 728 } else { 729 qmn = q1 + 0.5*Qstp; 730 qmx = q2 - 0.5*Qstp; 731 qst = (qmx-qmn)/floor((qmx-qmn)/Qstp+0.5); 732 qmx+= 1e-6*qst; 733 }; 734 for (qqq=qmn; qqq<=qmx; qqq+=qst) { 735 tst = 2.0*Pic/ddmax(floor(4.0*Pic*qqq/Qstp+0.5),12.0); 736 tmn = tst* floor( Tmn/tst-1.0); 737 tmx = -tst*(floor(-Tmx/tst-1.0)-1e-6); 738 if (tmx-tmn>=2.0*Pic) tmx = tmn + 2.0*Pic + 1e-6*tst; 739 for (ttt=tmn; ttt<=tmx; ttt+=tst) { 740 Qxx = qqq*cos(ttt)-Qx; 741 Qyy = qqq*sin(ttt)-Qy; 742 Ni = fabs(N11*Qxx+N12*Qyy); 743 Nj = fabs(N21*Qxx+N22*Qyy); 744// if (j==30 && i==30 && qi==qq1 && qqq==qmn && ttt==tmn && LAM>4.08 && LAM<4.30) fprintf(stdout," %g %g %g %g %g %g %g %g %g %g %g \n",QQQ,Qxx,Qyy,Ni,Nj,qmn,qmx,qst,tmn,tmx,tst); 745 if (Ni<1.0 && Nj<1.0) { 746// kkk++; 747 AR2 = (1.0-Ni)*(1.0-Nj); 748 FAK = prIn * Sthckn * OMG; 749 NNN = AR2*TSdN(i,j,s,k); 750 Int = AR2*TSdp(i,j,s,k) / FAK; 751 Err = AR2*TSdp2(i,j,s,k)/(FAK*FAK); 752 Wt = AR2*AREA; 753 754 calibN(qi,k) += NNN; 755 calibI(qi,k) += Int; 756 calibIe(qi,k)+= Err; 757 calibWt(qi,k)+= Wt; 758 calibN(qi,3) += NNN; 759 calibI(qi,3) += Int; 760 calibIe(qi,3)+= Err; 761 calibWt(qi,3)+= Wt; 762 763 si = floor(s3cl*20.0/Ntc); 764// if (si>=20 || si<0) fprintf(stdout,"Error in time"); 765 sectnN(qi,si) += NNN; 766 sectnI(qi,si) += Int; 767 sectnIe(qi,si)+= Err; 768 sectnWt(qi,si)+= Wt; 769 770 }; 771 }; 772 }; 773 }; // here check if primary intensity was enough 774 775// if (j==30 && i==30 && LAM>4.08 && LAM<4.30) fprintf(stdout,"%g %g %g %i \n",LAM,OMG,QQQ,kkk); 776 777 }; 778 }; 779 }; 780 }; 781 }; 782 kfl = 1; 783 }; 784 }; 785 786// fprintf(stdout,"point4"); 787 788 for (qi=0; qi<Nqc; qi++) { 789 for (si=0; si<4; si++) { 790 if (calibI(qi,si)>0.0 && calibWt(qi,si)>0.0) { 791// if (qi==110) fprintf(stdout,"%g\n",calibWt(qi,si)); 792 calibI(qi,si) /= calibWt(qi,si); 793 calibIe(qi,si)/= calibWt(qi,si)*calibWt(qi,si); 794 }; 795 }; 796 for (si=0; si<20; si++) { 797 if (sectnI(qi,si)>0.0 && sectnWt(qi,si)>0.0) { 798 sectnI(qi,si) /= sectnWt(qi,si); 799 sectnIe(qi,si)/= sectnWt(qi,si)*sectnWt(qi,si); 800 }; 801 }; 802 }; 803 804// fprintf(stdout,"point5"); 805 806 i=0; 807 while(i<99 && fname[i]>0) {filename[i]=fname[i]; i++;} 808 i=iimin(i,87); 809 j=i; 810 filename[j] = '_'; j++; 811 filename[j] = 'c'; j++; 812 filename[j] = 'p'; j++; 813 filename[j] = 's'; j++; 814 filename[j] = '_'; j++; 815 filename[j] = 'a'; j++; 816 filename[j] = 'l'; j++; 817 filename[j] = 'l'; j++; 818 filename[j] = '.'; j++; 819 filename[j] = 'd'; j++; 820 filename[j] = 'a'; j++; 821 filename[j] = 't'; j++; 822 filename[j] = 0; 823 824 DETECTOR_OUT_1D( 825 "TOFSANSdet.comp", 826 "log(Q) [AA^(-1)]", 827 "I(Q) [cps]", 828 "log(Q) [AA^(-1)]", qminl, qmaxl, Nq, 829 &simplN[0], &simplI[0], &simplIe[0], filename 830 ); 831 832 for (si=0; si<4; si++) { 833 834 j=i+1; 835 filename[j] = 'c'; j++; 836 filename[j] = 'a'; j++; 837 filename[j] = 'l'; j++; 838 filename[j] = 'i'; j++; 839 filename[j] = 'b'; j++; 840 filename[j] = '_'; j++; 841 filename[j] = 49+si; 842 843 for (qi=0; qi<Nqc; qi++) { 844// simplQ[qi]= calibQ(qi,si); 845 simplN[qi]= calibN(qi,si); 846 simplI[qi]= calibI(qi,si); 847 simplIe[qi]=calibIe(qi,si); 848 }; 849 850 DETECTOR_OUT_1D( 851 "TOFSANSdet.comp", 852 "log(Q) [AA^(-1)]", 853 "I(Q) [cm-1]", 854 "log(Q) [AA^(-1)]", qminl, qmaxl, Nq, 855 &simplN[0], &simplI[0], &simplIe[0], filename 856 ); 857 }; 858 859 860 for (si=0; si<20; si++) { 861 862 j=i+1; 863 filename[j] = 's'; j++; 864 filename[j] = 'e'; j++; 865 filename[j] = 'c'; j++; 866 filename[j] = 't'; j++; 867 filename[j] = '_'; j++; 868 filename[j] = 48+(si/10); j++; 869 filename[j] = 48+(si-10*(si/10)); 870 871 872 for (qi=0; qi<Nqc; qi++) { 873// simplQ[qi]= sectnQ(qi,si); 874 simplN[qi]= sectnN(qi,si); 875 simplI[qi]= sectnI(qi,si); 876 simplIe[qi]=sectnIe(qi,si); 877 }; 878 879 DETECTOR_OUT_1D( 880 "TOFSANSdet.comp", 881 "log(Q) [AA^(-1)]", 882 "I(Q) [cm-1]", 883 "log(Q) [AA^(-1)]", qminl, qmaxl, Nq, 884 &simplN[0], &simplI[0], &simplIe[0], filename 885 ); 886 }; 887 %} 888 889 890 891MCDISPLAY 892%{ 893 double ds1c,ds2c,ds3c; 894 double xw1c,xw2c,xw3c; 895 double yh1c,yh2c,yh3c; 896 double hl1c,hl2c,hl3c; 897 double vx3c; 898 899 ds1c = ddmax(ds1,0.0); 900 ds2c = ddmax(ds2,0.0); 901 ds3c = ddmax(ds3,0.0); 902 if (ds2c>=ds3c) ds2c = 0.0; 903 if (ds1c>=ds2c) ds1c = 0.0; 904 xw1c = ddmax(xw1,0.0); 905 yh1c = ddmax(yh1,0.0); 906 xw2c = ddmax(xw2,0.0); 907 yh2c = ddmax(yh2,0.0); 908 xw3c = ddmax(xw3,0.0); 909 yh3c = ddmax(yh3,0.0); 910 hl1c = hl1; 911 hl2c = hl2; 912 hl3c = hl3; 913 if (ds3c==0.0) {xw3c=0.0; yh3c=0.0; hl3c=0.0;}; 914 if (ds2c==0.0) {xw2c=0.0; yh2c=0.0; hl2c=0.0;}; 915 if (ds1c==0.0) {xw1c=0.0; yh1c=0.0; hl1c=0.0;}; 916 if (hl1c<=0.0) {hl1c = 0.0;} else {hl1c=ddmax(hl1c,ddmax(xw1c/Nx,yh1c/Ny));}; 917 if (hl2c<=0.0) {hl2c = 0.0;} else {hl2c=ddmax(hl2c,ddmax(xw2c/Nx,yh2c/Ny));}; 918 hl3c=ddmax(hl3c,ddmax(xw3c/Nx,yh3c/Ny)*3.0); /* leave some space for other data */ 919 vx3c = ddmin(ddmax(vx3,1.0),15.0); 920 921 922 923 if (ds1c>0.0 && xw1c>0.0 && yh1c>0.0) { 924 multiline(5, -0.5*xw1c, -0.5*yh1c, ds1c, 925 0.5*xw1c, -0.5*yh1c, ds1c, 926 0.5*xw1c, 0.5*yh1c, ds1c, 927 -0.5*xw1c, 0.5*yh1c, ds1c, 928 -0.5*xw1c, -0.5*yh1c, ds1c); 929 multiline(5, -0.5*hl1c, -0.5*hl1c, ds1c, 930 0.5*hl1c, -0.5*hl1c, ds1c, 931 0.5*hl1c, 0.5*hl1c, ds1c, 932 -0.5*hl1c, 0.5*hl1c, ds1c, 933 -0.5*hl1c, -0.5*hl1c, ds1c); }; 934 935 if (ds2c>0.0 && xw2c>0.0 && yh2c>0.0) { 936 multiline(5, -0.5*xw2c, -0.5*yh2c, ds2c, 937 0.5*xw2c, -0.5*yh2c, ds2c, 938 0.5*xw2c, 0.5*yh2c, ds2c, 939 -0.5*xw2c, 0.5*yh2c, ds2c, 940 -0.5*xw2c, -0.5*yh2c, ds2c); 941 multiline(5, -0.5*hl2c, -0.5*hl2c, ds2c, 942 0.5*hl2c, -0.5*hl2c, ds2c, 943 0.5*hl2c, 0.5*hl2c, ds2c, 944 -0.5*hl2c, 0.5*hl2c, ds2c, 945 -0.5*hl2c, -0.5*hl2c, ds2c); }; 946 947 if (ds3c>0.0 && xw3c>0.0 && yh3c>0.0) { 948 multiline(5, -0.5*xw3c, -0.5*yh3c, ds3c, 949 0.5*xw3c, -0.5*yh3c, ds3c, 950 0.5*xw3c, 0.5*yh3c, ds3c, 951 -0.5*xw3c, 0.5*yh3c, ds3c, 952 -0.5*xw3c, -0.5*yh3c, ds3c); 953 multiline(5, -0.5*hl3c, (0.5-vx3c)*hl3c, ds3c, 954 0.5*hl3c, (0.5-vx3c)*hl3c, ds3c, 955 0.5*hl3c, 0.5 *hl3c, ds3c, 956 -0.5*hl3c, 0.5 *hl3c, ds3c, 957 -0.5*hl3c, (0.5-vx3c)*hl3c, ds3c); }; 958%} 959 960 END 961