1/******************************************************************************* 2* 3* McStas, neutron ray-tracing package 4* Copyright (C) 1997-2011, All rights reserved 5* Risoe National Laboratory, Roskilde, Denmark 6* Institut Laue Langevin, Grenoble, France 7* 8* Component: curved_mirror 9* 10* %I 11* Written by: Henrik Jacobsen 12* Date: April 2012 13* Origin: RNBI 14* 15* Single mirror plate that is curved and fits into an elliptic guide. 16* 17* %D 18* 19* 20* %P 21* INPUT PARAMETERS: 22* 23* focus_s: [m] first focal point of ellipse in ABSOLUTE COORDINATES 24* focus_e: [m] second focal point of ellipse in ABSOLUTE COORDINATES 25* mirror_start: [m] start of mirror in ABSOLUTE COORDINATES 26* guide_start: [m] start of guide in ABSOLUTE COORDINATES 27* yheight: [m] the height of the mirror when not in the guide 28* smallaxis: [m] the smallaxis of the guide 29* length: [m] the length of the mirror 30* R0=0.99: [] low angle reflectivity 31* Qc: [AA-1] critical scattering vector 32* alpha: [AA] slope of reflectivity 33* m: [] m-value of material 34* W: [AA-1] width of supermirror cutoff 35* transmit: [0/1] 0: non reflected neutrons are absorbed. 1: non reflected neutrons pass through 36* substrate_thickness: [m] thickness of substrate (for absorption) 37* coating_thickness: [m] thickness of coating (for absorption) 38* theta_1 [deg] angle of mirror wrt propagation direction at start of mirror 39* theta_2 [deg] angle of mirror wrt propagation direction at center of mirror 40* theta_3 [deg] angle of mirror wrt propagation direction at end of mirror 41* reflect: [q(Angs-1) R(0-1)] (str) Name of relfectivity file. Format 42* 43* %D 44* 45* %E 46*******************************************************************************/ 47 48DEFINE COMPONENT Mirror_Curved_Bispectral 49DEFINITION PARAMETERS (string reflect=0) 50SETTING PARAMETERS (focus_s,focus_e, mirror_start, guide_start, yheight, smallaxis, length ,m, transmit=0, substrate_thickness=0.0005, coating_thickness=10e-6, theta_1=1.25, theta_2=1.25, theta_3=1.25) 51OUTPUT PARAMETERS (pTable) 52 53SHARE 54%{ 55%include "read_table-lib" 56%} 57 58DECLARE 59%{ 60t_Table pTable; 61 62double f; //half of distance between focal points 63double asquared; 64double a; //half of ellipse length 65double b; //half of ellipse width 66 67double xprime; //x in coordinates with center of ellipse at (xprime,zprime)=(0,0) 68double ymirror; //height of the mirror 69 70 71//Defining the mirror 72double a1; 73double b1; 74double c1; 75 76//solving the time the neutron hits the sample 77double A, B, C, D, E, P, Q, R, U, V, I, J, K; 78 79//finding rotation of mirror 80double alpha1, beta1, gamma1; 81double theta_m; 82double sin_theta_m, cos_theta_m; 83 84double tan_theta_1; 85double tan_theta_2; 86double tan_theta_3; 87 88 89double v_n; //speed of neutron perpendicular to surface 90 91double Ref; //reflectivity 92 93double dt; 94double q; 95 int intersect; 96 97double discriminant; 98 99 100 101 102double dt_2; 103double dt_3; 104int prop_case; 105double x_2; 106double y_2; 107double z_2; 108double t_2; 109double x_3; 110double y_3; 111double z_3; 112double t_3; 113 114int x_hit; 115int x_hit_2; 116int x_hit_3; 117double xprime_2; 118double ymirror_2; 119double xprime_3; 120double ymirror_3; 121int intersect_2; 122int intersect_3; 123 124 125 126%} 127 128INITIALIZE 129%{ 130 if (reflect && strlen(reflect) && strcmp(reflect,"NULL") && strcmp(reflect,"0")) { 131 if (Table_Read(&pTable, reflect, 1) <= 0) /* read 1st block data from file into pTable */ 132 exit(fprintf(stderr,"Mirror: %s: can not read file %s\n", NAME_CURRENT_COMP, reflect)); 133 } 134%} 135 136TRACE 137%{ 138intersect=0; 139x_hit=0; 140x_hit_2=0; 141x_hit_3=0; 142intersect_2=0; 143intersect_3=0; 144prop_case=0; 145 146//printf("\n\n\n"); 147 double old_x = x, old_y = y, old_z = z, old_t=t, old_vx=vx, old_vz=vz, old_vy=vy; 148 149// printf("x=%f, y=%f, z=%f, vx=%f, vy=%f, vz=%f\n",x,y,z,vx,vy,vz); 150 151// Check if neutron hits mirror. First find which z,x coordinates it hits. 152 153//mirror is defined by z(x)=a1x^3+b1x^2+c1x+d1, with dz/dx|x=-length/2=tan(theta_1), dz/dx|x=0=tan(theta_2), dz/dx|x=length/2=tan(theta3), z(0)=0. (d1=0) 154 155tan_theta_1=tan(theta_1*DEG2RAD); 156tan_theta_2=tan(theta_2*DEG2RAD); 157tan_theta_3=tan(theta_3*DEG2RAD); 158 159 160a1=2.0/3.0*(tan_theta_1+tan_theta_3-2.0*tan_theta_2)/(length*length); 161b1=(tan_theta_3-tan_theta_1)/(2.0*length); 162c1=tan_theta_2; 163 164 165//neutron trajectory is defined by x=x0+vx*t, z=z0+vz*t. setting z=a1*x^3+b1*x^2+c1*x gives the equation A*t^3+B*t^2+C*t+D=0, with 166A=a1*vx*vx*vx; 167B=3.0*a1*x*vx*vx+b1*vx*vx; 168C=3.0*a1*x*x*vx+2.0*b1*x*vx+c1*vx-vz; 169D=a1*x*x*x+b1*x*x+c1*x-z; 170 171//printf("a1=%f,b1=%f,c1=%f",a1,b1,c1); 172 173//this equation must now be solved for t; 174 175if (A!=0){ 176P=1/3.0*(3.0*C/A-B*B/(A*A)); 177Q=1/27.0*(2.0*B*B*B/(A*A*A)-9.0*B*C/(A*A)+27.0*D/A); 178 179E=P*P*P/27.0+Q*Q/4.0; 180 181// printf("A=%f, B=%f, C=%f, D=%f, 1e6P=%f, 1e6Q=%f, 1e6E=%f\n", A, B, C, D, 1e6*P, 1e6*Q, 1e6*E); 182 183prop_case=0; 184if (E>=0){ 185 186U=cbrt(-Q/2.0+sqrt(E)); 187V=cbrt(-Q/2.0-sqrt(E)); 188 189I=U+V-B/(3.0*A); 190dt=I; 191dt_2=I; 192dt_3=I; 193// printf("I=%f\n",I); 194 195// J=-(U+V)/2+1i*(U-V)*sqrt(3)/2-B/(3*A) //complex solution 196// K=-(U+V)/2-1i*(U-V)*sqrt(3)/2-B/(3*A) //complex solution 197}else{ 198 R=acos(-Q/(2.0*sqrt(-P*P*P/27.0))); 199 200// printf("R=%f\n",R); 201 202 203 I=2.0*sqrt(fabs(P)/3.0)*cos(R/3.0)-B/A/3.0; 204 J=-2.0*sqrt(fabs(P)/3.0)*cos(R/3.0 + 3.1415926535/3.0)-B/A/3.0; 205 K=-2.0*sqrt(fabs(P)/3.0)*cos(R/3.0 - 3.1415926535/3.0)-B/A/3.0; 206 207// printf("2.0*sqrt(abs(P)/3.0)=%f", 2.0*sqrt(abs(P)/3.0)); 208// printf("cos(R/3.0)=%f, cos(R/3.0 + 3.1415926535/3.0)=%f, cos(R/3.0 - 3.1415926535/3.0)=%f, -B/A/3.0=%f\n", cos(R/3.0), cos(R/3.0 + 3.1415926535/3.0), cos(R/3.0 - 3.1415926535/3.0), -B/A/3.0); 209 210// printf("I=%f, J=%f, K=%f, \n",I, J, K); 211// printf("P=%f, R=%f, A=%f, B=%f, \n",P, R, A, B); 212 213 214 215// Three solutions. Find the smallest positive of these. 216//there are problems with the solutions.... 217 if (I<=0){ 218 if (J<=0 && K<=0){dt=-1.0;} //if all three are negative, dt<0 and nothing happens 219 if (J<=0 && K>0){dt=K;} //if only K>0, dt=K 220 if (J>0 && K<=0){dt=J;} //if only J>0, dt=J 221 222 if (J>0 && K>0){ //if both J>0 and K>0, compare 223 if (J>=K){dt=K; prop_case=1; dt_2=J;}else{dt=J; dt_2=I; prop_case=2;} } //dt is the smallest value 224 }else{ //end if (I<=0) 225 if (J<=0 && K<=0){dt=I;} //if only I>0, dt=I; 226 227 if (J<=0 && K>0){ //if both I>0 and K>0, compare 228 if (K>=I){dt=I; dt_2=K; prop_case=3;}else{dt=K; dt_2=I; prop_case=4;} } //dt is the smallest value 229 230 if (J>0 && K<=0){ //if both I>0 and J>0, compare 231 if (J>=I){dt=I; dt_2=J; prop_case=5;}else{dt=J; dt_2=I; prop_case=6;} } //dt is the smallest value 232 233 if (J>0 && K>0){ //if all three>0, compare 234 if (J>=K){ //either K or I is smallest 235 if (K>=I){dt=I; if(J>=K){ dt_2=K; dt_3=J; prop_case=9;}else{dt_2=K; dt_3=J; prop_case=15;}}else{dt=K; if (J>=I){dt_2=I; dt_3=J; prop_case=10;}else{dt_2=J; dt_3=I; prop_case=11;} } //if K is smallest, compare it to I 236 }else{ 237 if (J>=I){dt=I; if (K>J){dt_2=J; dt_3=K; prop_case=12;}else{dt_2=J; dt_3=J; prop_case=16;}}else{dt=J; if (K>I){dt_2=I; dt_3=K; prop_case=13;}else{{dt_2=K; dt_3=I; prop_case=14;}} }} //else compare J to I 238 } //end if(J>0 && K>0) 239 240 } //end }else{ for if(I<=0) 241 242 243 244} // end }else{ for if (E>=0) 245 246 247}else{ //end if (A!=0) 248if (B!=0){ 249 250discriminant=C*C-4*B*D; 251 252if (discriminant<0){dt=-1.0;}else{ //only complex solutions: set dt<0 to avoid interaction 253I=(-C-sqrt(discriminant))/(2.0*B); 254J=(-C+sqrt(discriminant))/(2.0*B); 255 256if (I<=0 && J<=0){dt = -1.0;} //both times are negative. 257if (I<=0 && J>0 ){dt = J;} //set dt to only positive value. 258if (I>0 && J<=0){dt = I;} //set dt to only positive value. 259if (I>0 && J>0 ){if (I>J) {dt=J; dt_2=I; prop_case=7;}else{dt=I; dt_2=J; prop_case=8;} } //set dt to smallest positive value 260 261} //end if (discriminant<0){}else{ 262}else{ //end if (B!)=0 263if (C!=0) { dt = -D/C;}else{ 264 printf("warning: A=B=C=0. Neutron is ignored\n"); } 265} //end if(B!=0){}else{ 266} //end if (A!=0){}else{ 267//now intersection time has been found. 268 269if (dt>0) { //if time is positive, propagate neutron to where it hits mirror. This is done without gravity. 270// printf("before anything: x=%f,y=%f,z=%f,vx=%f,vy=%f,vz=%f, dt=%f\n",x,y,z,vx,vy,vz,dt); 271 272 x += vx*dt; 273 y += vy*dt; 274 z += vz*dt; 275 t += dt; 276 277 278x_hit=(x >=-length/2 && x<=length/2); 279 280 281if (prop_case==0){ 282x_2=x; 283y_2=y; 284z_2=z; 285t_2=t; 286x_3=x; 287y_3=y; 288z_3=z; 289t_3=t; 290} 291 292if (prop_case>0) 293{ 294x_2=old_x+vx*dt_2; 295y_2=old_y+vy*dt_2; 296z_2=old_z+vz*dt_2; 297t_2=old_t+dt_2; 298x_hit_2=(x_2 >=-length/2 && x_2<=length/2); 299} 300 301if (prop_case>8) 302{ 303x_3=old_x+vx*dt_3; 304y_3=old_y+vy*dt_3; 305z_3=old_z+vz*dt_3; 306t_3=old_t+dt_3; 307x_hit_3=(x_3 >=-length/2 && x_3<=length/2); 308} 309 310//printf("x_hit=%d, x_hit_2=%d, x_hit_3=%d\n",x_hit, x_hit_2, x_hit_3); 311//printf("dt=%f, dt_2=%f, dt_3=%f\n",dt,dt_2,dt_3); 312// printf("x=%f,y=%f,z=%f,vx=%f,vy=%f,vz=%f\n",x,y,z,vx,vy,vz); 313 314// printf("x=%f, length/2=%f\n",x, length/2); 315 316 317if (x_hit || x_hit_2 || x_hit_3){ 318//if (x >=-length/2 && x<=length/2){ //check if neutron is within x limits of the mirror. If so, check if it is within y limits. 319 320 321//define the ellipse 322b=smallaxis/2; 323 324f=(focus_e-focus_s)*0.5; 325 326 asquared=f*f+b*b; 327 a=sqrt(asquared); 328 329xprime=-f-focus_s+mirror_start+length/2+x; //xprime is the x-coordinate in a coordinate system centered at the center of the ellipse 330 331//ymirror=b*sqrt(1-xprime*xprime/(f*f)); //following Kaspars convention, assuming f~=a (valid for most elliptic guides normally used) 332 333ymirror=b*sqrt(1-xprime*xprime/asquared); 334 335 336 337xprime_2=-f-focus_s+mirror_start+length/2+x_2; //xprime is the x-coordinate in a coordinate system centered at the center of the ellipse 338ymirror_2=b*sqrt(1-xprime_2*xprime_2/asquared); 339 340xprime_3=-f-focus_s+mirror_start+length/2+x_3; //xprime is the x-coordinate in a coordinate system centered at the center of the ellipse 341ymirror_3=b*sqrt(1-xprime_3*xprime_3/asquared); 342 343if (guide_start>mirror_start){ //If (part of the) mirror is outside the guide, the mirror can be extended 344if ( x<-length/2+guide_start-mirror_start) { 345ymirror=yheight/2; 346} 347 348if ( x_2<-length/2+guide_start-mirror_start) { 349ymirror_2=yheight/2; 350} 351 352if ( x_3<-length/2+guide_start-mirror_start) { 353ymirror_3=yheight/2; 354} 355 356 357 358 359} 360 361 362 363 364 365 366 367 368 369 370 371// printf("ymirror=%f, y=%f\n",ymirror, y); 372intersect = ( y>=-ymirror && y<=ymirror && x >=-length/2 && x<=length/2); 373 374if (prop_case>0) { 375intersect_2 = ( y_2>=-ymirror && y_2<=ymirror && x_2 >=-length/2 && x_2<=length/2); 376} 377if (prop_case>8){ 378intersect_3 = ( y_3>=-ymirror && y_3<=ymirror && x_3 >=-length/2 && x_3<=length/2); 379} 380 381//printf("y_2=%f, ymirror=%f\n",y_2,ymirror); 382 383//printf("\nintersect=%d, intersect_2=%d, intersect_3=%d, prop_case=%d\n",intersect, intersect_2, intersect_3, prop_case); 384 385//printf("x=%f,y=%f,z=%f,t=%f\n",x,y,z,t); 386//printf("x_2=%f,y_2=%f,z_2=%f,t_2=%f\n",x_2,y_2,z_2,t_2); 387//printf("x_3=%f,y_3=%f,z_3=%f,t_3=%f\n",x_3,y_3,z_3,t_3); 388 389if (!intersect){ 390if (!intersect_2){ 391intersect=intersect_3; 392x=x_3; 393y=y_3; 394z=z_3; 395t=t_3; 396}else{ 397intersect=intersect_2; 398x=x_2; 399y=y_2; 400z=z_2; 401t=t_2; 402} 403} 404 405//printf("intersect=%d, intersect_2=%d, intersect_3=%d, prop_case=%d\n\n",intersect, intersect_2, intersect_3, prop_case); 406//printf("x=%f,y=%f,z=%f,t=%f\n",x,y,z,t); 407 408//printf("z=%f, zcalc=%f\n",z,a1*x*x*x+b1*x*x+c1*x); 409//printf("z=%f, zcalc=%f\n",z_2,a1*x_2*x_2*x_2+b1*x_2*x_2+c1*x_2); 410//printf("z=%f, zcalc=%f\n",z_3,a1*x_3*x_3*x_3+b1*x_3*x_3+c1*x_3); 411 412 if (intersect) { //if neutron is within ylimits of the mirror handle reflection/transmission 413 414//first find the angle of the mirror. It is given by theta(x)=alpha*x^2+beta*x+gamma1, with theta(-l/2)=theta1, theta(0)=theta2, theta(l/2)=theta3 415 416alpha1=2*(theta_1+theta_3-2*theta_2)/(length*length); 417beta1=(theta_3-theta_1)/length; 418gamma1=theta_2; 419 420theta_m=alpha1*x*x+beta1*x+gamma1; // angle of mirror. 421 422//The vector normal to the mirror is e_n= sin(theta)*e_x-cos(theta)*e_z 423 424//find amplitude of v in direction of e_n: 425 426sin_theta_m=sin(theta_m*DEG2RAD); 427cos_theta_m=cos(theta_m*DEG2RAD); 428 429v_n=sin_theta_m*vx-cos_theta_m*vz; 430 431 432q=fabs(2.0*v_n*V2Q); 433 434double R0=0.99; 435double Qc=0.0217; 436double m_value=m*0.9853+0.1978; 437double W=-0.0002*m_value+0.0022; 438double alpha=0.1204*m_value+5.0944; 439double beta=-7.6251*m_value+68.1137; 440 441if (m_value<=3) 442{alpha=m_value; 443beta=0;} 444 445 446 447 448 449 /* Reflectivity (see component Guide). */ 450 if(m == 0) 451 ABSORB; 452 if (reflect && strlen(reflect) && strcmp(reflect,"NULL") && strcmp(reflect,"0")) 453 Ref=Table_Value(pTable, q, 1); 454 else { 455 Ref = R0; 456 if(q > Qc) 457 { 458 double arg = (q-m_value*Qc)/W; 459 if(arg < 10) 460 Ref *= .5*(1-tanh(arg))*(1-alpha*(q-Qc)+beta*(q-Qc)*(q-Qc)); //matches data from Swiss Neutronics 461 else Ref=0; 462 } 463 } 464 if (Ref < 0) Ref=0; 465 else if (Ref > 1) Ref=1; 466 467 468//Now comes actual reflection/transmission 469 if (!transmit) { //all neutrons are reflected 470 if (!Ref) ABSORB; 471 p *= Ref; 472 473//handle reflection: change v_n -->-v_n 474 475vx=old_vx*(cos_theta_m*cos_theta_m-sin_theta_m*sin_theta_m)+old_vz*(2*cos_theta_m*sin_theta_m); 476vz=old_vx*(2*cos_theta_m*sin_theta_m)+old_vz*(sin_theta_m*sin_theta_m-cos_theta_m*cos_theta_m); 477 478// printf("theta_m=%f, sin_theta_m=%f, cos_theta_m=%f, v_n=%f, old_vx=%f, vx=%f, old_vz=%f, vz=%f\n\n", theta_m, sin_theta_m, cos_theta_m, v_n, old_vx, vx, old_vz, vz); 479 480 481 SCATTER; 482//printf("line 471.In mirror: x=%f,y=%f,z=%f,t=%f\n",x,y,z,t); 483//printf("In mirror: old_vx=%f,old_vy=%f,old_vz=%f,vx=%f,vy=%f,vz=%f,v_n=%f\n",old_vx,old_vy,old_vz,vx,vy,vz,v_n); 484 485 } else { //if neutrons can be transmitted 486 487 488 489//calculate absorption. 490// substrate 491double lambda=(2*PI/V2K)/sqrt(vx*vx + vy*vy + vz*vz); 492double sin_theta=lambda*q/(4*PI); 493 494//double substrate_path_length=substrate_thickness/sin_theta; 495//double coating_path_length=coating_thickness/sin_theta; 496 497double sin_theta_c=Qc/(4*PI); 498 499double theta_diff; 500double substrate_path_length; 501double coating_path_length; 502 503double remaining_length_through_mirror; 504 505int hit_back_mirror; 506 507if (v_n>0) { 508hit_back_mirror=1;} else{ 509hit_back_mirror=0;} 510 511remaining_length_through_mirror=length/2-x; 512 513 514if (sin_theta>sin_theta_c*lambda) { 515theta_diff=sqrt(sin_theta*sin_theta-sin_theta_c*sin_theta_c*lambda*lambda); 516coating_path_length=coating_thickness/theta_diff; 517substrate_path_length=substrate_thickness/theta_diff; 518 519 if (coating_path_length>remaining_length_through_mirror){ 520coating_path_length=remaining_length_through_mirror; 521substrate_path_length=0; 522} 523 524 if (substrate_path_length>remaining_length_through_mirror){ 525substrate_path_length=remaining_length_through_mirror; 526} 527 528 529 530 531 532 533 534 535 536 537 538 539} else{ 540 541if (hit_back_mirror==0){ //neutron comes from front of mirror 542substrate_path_length=0; 543coating_path_length=remaining_length_through_mirror; 544}else {//neutron comes from behind mirror 545 546substrate_path_length=remaining_length_through_mirror; 547coating_path_length=0; 548} 549} 550 551 552double mu_substrate=0.0318/lambda+0.0055*lambda-0.0050; //unit: cm^-1 553mu_substrate=mu_substrate*100; //unit: m^-1; 554 555//For nickel and titanium coating, the following formular is used: 556// mu = rho/m(atom)*sigma_a,thermal*lambda/lambda_thermal 557 558// lambda_thermal=1.798 Å 559 560// rho_nickel=8.908g/cm^3 561// m(atom)_nickel=58.6934*1.661*10^-27 kg 562// sigma_a,thermal_nickel=4.49*10^-28 m^2 563 564// rho_titanium=4.506g/cm^3 565// m(atom)_titanium=47.867*1.661*10^-27 kg 566// sigma_a,thermal_titanium=6.09*10^-28 m^2 567 568double Ni_coefficient=22.8180; 569double Ti_coefficient=19.1961; 570 571double mu_coating=(0.5*Ni_coefficient+0.5*Ti_coefficient)*lambda; //it is roughly 50% nickel and 50% titanium 572 573 574 575 // transmit when rand > R 576 if (Ref == 0 || rand01() >= Ref) { //transmit 577if (substrate_thickness>0){ p=p*exp(-mu_substrate*substrate_path_length-mu_coating*coating_path_length); //reduce weight of neutrons due to attenuation in the mirror 578//x+=(coating_path_length+substrate_path_length)-(coating_thickness+substrate_thickness)/sin_theta; 579//printf("xshift is %f \n",(coating_path_length+substrate_path_length)-(coating_thickness+substrate_thickness)/sin_theta); 580} 581// printf("line 380\n"); 582/* 583if (v_n>0) { 584printf("neutron is transmitted from back of mirror. %f\n",exp(-mu_substrate*substrate_path_length-mu_coating*coating_path_length)); 585}else{ 586printf("neutron is transmitted from front of mirror. %f\n",exp(-mu_substrate*substrate_path_length-mu_coating*coating_path_length)); 587} 588*/ 589} else {//neutron is reflected 590 if (v_n>0 && substrate_thickness>0) { //if neutron comes from behind the mirror 591// printf("neutron is reflected from back of mirror. %f\n",Ref*exp(-2*mu_substrate*substrate_path_length-2*mu_coating*coating_path_length)); 592 p=p*exp(-2*mu_substrate*substrate_path_length-2*mu_coating*coating_path_length);} //else{ //reduce weight of neutrons due to attenuation in the mirror 593 // printf("neutron is reflected from front of mirror. %f\n", Ref);} 594//handle reflection: change v_n -->-v_n 595vx=old_vx*(cos_theta_m*cos_theta_m-sin_theta_m*sin_theta_m)+old_vz*(2*cos_theta_m*sin_theta_m); 596vz=old_vx*(2*cos_theta_m*sin_theta_m)+old_vz*(sin_theta_m*sin_theta_m-cos_theta_m*cos_theta_m); 597// printf("line 388\n"); 598 599} 600 601// printf("theta_m=%f, sin_theta_m=%f, cos_theta_m=%f, v_n=%f, old_vx=%f, vx=%f, old_vz=%f, vz=%f\n\n", theta_m, sin_theta_m, cos_theta_m, v_n, old_vx, vx, old_vz, vz); 602 603//printf("vxvx+vzvz=%f, oldvxoldvx+oldvzoldvz=%f", vx*vx+vz*vz, old_vx*old_vx+old_vz*old_vz); 604 605 SCATTER; 606//printf("line 524.In mirror: x=%f,y=%f,z=%f,t=%f\n",x,y,z,t); 607//printf("old_vx=%f,old_vy=%f,old_vz=%f,vx=%f,vy=%f,vz=%f,v_n=%f\n",old_vx,old_vy,old_vz,vx,vy,vz,v_n); 608//after transmission or reflection 609 } //end } else { after if (!transmit) { 610 } 611 612 613 614 615 616} // end if (x >=-length/2 && x<=length/2) 617 618// printf("intersect=%d\n",intersect); 619 620 if (!intersect) { 621 /* No intersection: restore neutron position. */ 622 x = old_x; 623 y = old_y; 624 625 z = old_z; 626 t = old_t; 627// printf("line 409\n"); 628 629 } 630 631 632} //end if (dt>0) 633 634 635%} 636 637MCDISPLAY 638%{ 639 640 641 642/* 643if (xcenter==0){ 644 645xstart=0; 646xend=length; 647xprime_start=-a+mirror_start+xstart; 648ystart=b*sqrt(1-xprime_start*xprime_start/asquared); 649 650xprime_end=-a+mirror_start+xend; 651yend=b*sqrt(1-xprime_end*xprime_end/asquared); 652 653} 654*/ 655 656/* 657if (xcenter==1){ 658xstart=-length/2; 659xend=length/2; 660 661 662xprime_start=-a+mirror_start+xstart+length/2; 663ystart=b*sqrt(1-xprime_start*xprime_start/asquared); 664 665xprime_end=-a+mirror_start+xend+length/2; 666yend=b*sqrt(1-xprime_end*xprime_end/asquared); 667} 668 669line(xstart,-ystart,0,xstart,ystart,0); 670line(xend,-yend,0,xend,yend,0); 671line(xstart,-ystart,0,xend,-yend,0); 672line(xstart,ystart,0,xend,yend,0); 673*/ 674 675double xstart; 676double xend; 677double xprime_start; 678double ystart; 679 680double xprime_end; 681double yend; 682 683double focus_2; 684double focus_1; 685double b; 686double f; 687double asquared; 688double a; 689 690 691int n_lines; 692int j=0; 693double xprimepos[51]; 694double ypos[51]; 695double x_plot[51]; 696double zpos[51]; 697double xstep; 698 699 700 701focus_2=focus_e-mirror_start; //focus in local coordinates 702focus_1=focus_s-mirror_start; 703 704b=smallaxis/2; 705 706f=(focus_2-focus_1)*0.5; 707 asquared=f*f+b*b; 708 a=sqrt(asquared); 709 710 711xstart=-length/2; 712xend=length/2; 713 714n_lines=50; 715 716xstep=length/n_lines; 717 718 719 720 721 722 723 724double a1, b1, c1; 725double tan_theta_1; 726double tan_theta_2; 727double tan_theta_3; 728 729 730 731//mirror is defined by z(x)=a1x^3+b1x^2+c1x+d1, with dz/dx|x=-length/2=tan(theta_1), dz/dx|x=0=tan(theta_2), dz/dx|x=length/2=tan(theta3), z(0)=0. (d1=0) 732 733tan_theta_1=tan(theta_1*DEG2RAD); 734tan_theta_2=tan(theta_2*DEG2RAD); 735tan_theta_3=tan(theta_3*DEG2RAD); 736 737 738a1=2.0/3.0*(tan_theta_1+tan_theta_3-2.0*tan_theta_2)/(length*length); 739b1=(tan_theta_3-tan_theta_1)/(2.0*length); 740c1=tan_theta_2; 741 742 743 744for (j=0; j<n_lines+1; j++) 745{ 746xprimepos[j]=-f-focus_s+mirror_start+length/2+xstart+xstep*j; 747 748ypos[j]=b*sqrt(1-xprimepos[j]*xprimepos[j]/asquared); //correct 749 750if (guide_start>mirror_start){ 751if ( xstart+xstep*j<-length/2+guide_start-mirror_start) { 752ypos[j]=yheight/2; 753} 754} 755 756 757 758// ypos[j]=b*sqrt(1-xprimepos[j]*xprimepos[j]/(f*f)); //following convention in Kaspar's elliptic guide.. 759// printf("xprimepos[j]=%f,f*f=%f, ypos[j]=%f\n",xprimepos[j],f*f,ypos[j]); 760 761x_plot[j]=xstart+xstep*j; 762zpos[j]=a1*x_plot[j]*x_plot[j]*x_plot[j]+b1*x_plot[j]*x_plot[j]+c1*x_plot[j]; 763} 764 765for (j=0; j<n_lines; j++) 766{ 767line(x_plot[j], -ypos[j], zpos[j], x_plot[j+1], -ypos[j+1],zpos[j+1]); 768line(x_plot[j], ypos[j], zpos[j], x_plot[j+1], ypos[j+1],zpos[j+1]); 769} 770 771 772line(x_plot[0],-ypos[0],zpos[0],x_plot[0],ypos[0],zpos[0]); 773line(x_plot[50],-ypos[50],zpos[50],x_plot[50],ypos[50],zpos[50]); 774 775 776 777 778//printf("ypos0=%f xpos0=%f ypos50=%f xpos50=%f",ypos[0], x_plot[0], ypos[50], x_plot[50]); 779 780/* double xmax, xmin, ymax, ymin; 781 782 783 if (center == 0) { 784 xmax= x1; xmin=0; 785 ymax= yheight; ymin=0; 786 } else { 787 xmax= x1/2; xmin=-xmax; 788 ymax= yheight/2; ymin=-ymax; 789 } 790 multiline(5, (double)xmin, (double)ymin, 0.0, 791 (double)xmax, (double)ymin, 0.0, 792 (double)xmax, (double)ymax, 0.0, 793 (double)xmin, (double)ymax, 0.0, 794 (double)xmin, (double)ymin, 0.0); 795*/ 796%} 797END 798