/******************************************************************************* * * McStas, the neutron ray-tracing package: Guide_tapering.comp * * Component: Guide_tapering * * %I * Written by: Uwe Filges * Date: 22/10/2003 * Version: $Revision: 1.6 $ * Origin: McStas 1.6 * * Models a rectangular tapered guide (many shapes) * * %D * Models a rectangular guide tube centered on the Z axis. The entrance lies * in the X-Y plane. * The guide may be tapered. * * The component includes a feature to read in self-defined functions for * guide tapering. Under the parameter 'option' the KEYWORD 'file=' offers * the possibility to read in parameters from an ASC-file. The file structure * is shown below in the example. It is important to know that the first * 3 lines will be interpreted as comments. * Afterwards the dimension of each guide segment must be defined. The * length of each segment is constant l(i)=l/segments . The number of * segments is defined through the number of lines minus the first 3 * lines taken from the Input-File. * * Example Input-File: * * c Guide_tapering.comp * c i = 0 - 199 segments * c h1(i) h2(i) w1(i) w2(i) * 0.120000 0.119850 0.020000 0.020000 * 0.119850 0.119700 0.020000 0.020000 * 0.119700 0.119550 0.020000 0.020000 * 0.119550 0.119400 0.020000 0.020000 * 0.119400 0.119250 0.020000 0.020000 * 0.119250 0.119100 0.020000 0.020000 * ... * * Example1:Guide_tapering(w1=0.1, h1=0.18, linw=0.1, loutw=0.1, * linh=0.1, louth=0.1, l=1.5, option="elliptical", R0=0.99, * Qcx=0.021, Qcy=0.021, alphax=6.07, alphay=6.07, W=0.003, * mx=1, my=1, segno=800) * * Example2:Guide_tapering(w1=0, h1=0, linw=0, loutw=0, * linh=0, louth=0, l=1.5, option="file=ownfunction.txt", R0=0.99, * Qcx=0.021, Qcy=0.021, alphax=6.07, alphay=6.07, W=0.003, * mx=1, my=1) * * %BUGS * This component does not work with gravitation on. Use component Guide_gravity then. * * %P * INPUT PARAMETERS: * * w1: [m] Width at the guide entry * h1: [m] Height at the guide entry * linw [m] distance of 1. focal point and real guide entry - * left and right horizontal mirrors * loutw [m] distance of real guide exit and 2nd focal point - * left and right horizontal mirrors * l: [m] length of guide * linh: [m] distance of 1st focal point and real guide entry - * top and bottom vertical mirrors * louth: [m] distance of real guide exit and 2nd focal point - * top and bottom vertical mirrors * option [str] define the input function for the curve of the guide walls * options are: * "elliptical" - define elliptical function of guide walls * "parabolical" - define parabolical function of guide walls * "file=[filename]" - read in ASC-file with arbitrary * definition for the curve of the guide walls * R0: [1] Low-angle reflectivity * Qcx: [AA-1] Critical scattering vector for left and right vertical * mirrors in each channel * Qcy: [AA-1] Critical scattering vector for top and bottom mirrors * alphax: [AA] Slope of reflectivity for left and right vertical * mirrors in each channel * alphay: [AA] Slope of reflectivity for top and bottom mirrors * mx: [1] m-value of material for left and right vertical mirrors * in each channel. Zero means completely absorbing. * my: [1] m-value of material for top and bottom mirrors. Zero * means completely absorbing. * W: [AA-1] Width of supermirror cut-off for all mirrors * segno: [1] number of segments (z-axis) for cutting the tube * * %End * ******************************************************************************/ DEFINE COMPONENT Guide_tapering DEFINITION PARAMETERS (option=0,segno=800) SETTING PARAMETERS (w1=0,h1=0,l,linw=0,loutw=0,linh=0,louth=0, R0=0.99, Qcx=0.021,Qcy=0.021, alphax=6.07, alphay=6.07, W=0.003, mx=1, my=1) OUTPUT PARAMETERS(w1c,w2c,ww,hh,whalf,hhalf,lwhalf,lhhalf,h1_in,h2_out,w1_in,w2_out,l_seg) STATE PARAMETERS (x,y,z,vx,vy,vz,t,sx,sy,p) DECLARE %{ double w1c[segno]; double w2c[segno]; double ww[segno], hh[segno]; double whalf[segno], hhalf[segno]; double lwhalf[segno], lhhalf[segno]; double h1_in[segno+1], h2_out[segno+1], w1_in[segno+1], w2_out[segno+1]; double l_seg, h12, h2, w12, w2, a_ell_q, b_ell_q, lbw, lbh; double mxi ,u1 ,u2 ,div1, p2_para; double test,Div1; int i,ii,seg; char *fu; char *pos; char file_name[1024]; char *ep; FILE *num; %} INITIALIZE %{ struct para { char st[128]; } segment[800]; if (W <=0) { fprintf(stderr,"Component: %s (Guide_tapering) W must \n", NAME_CURRENT_COMP); fprintf(stderr," be positive\n"); exit(-1); } if (l <= 0) { fprintf(stderr,"Component: %s (Guide_tapering) real guide length \n", NAME_CURRENT_COMP); fprintf(stderr," is <= ZERO ! \n"); exit(-1); } if (mcgravitation) fprintf(stderr,"WARNING: Guide_tapering: %s: " "This component produces wrong results with gravitation !\n" "Use Guide_gravity.\n", NAME_CURRENT_COMP); seg=segno; l_seg=l/seg; h12 = h1/2; if (option != NULL) { fu = (char*)malloc(sizeof(char)*(strlen(option)+1)); strcpy(fu,option); } else { exit(-1); } if (!strcmp(fu,"elliptical")) { /* calculate parameter b of elliptical equestion - vertical mirrors */ /* (l+linh+louth) -> distance between focal points */ /* printf("A1 \n"); */ lbh = l + linh + louth; if (linh == 0 && louth == 0 ) { /* plane mirrors (vertical) */ b_ell_q = 0; h2 = h1; } else { /* elliptical mirrors */ u1 = sqrt((linh*linh)+(h12*h12)); u2 = sqrt((h12*h12) + ((l+louth)*(l+louth))); a_ell_q = ((u1 + u2)/2)*((u1 + u2)/2); b_ell_q = a_ell_q - ((lbh/2)*(lbh/2)); /* calculate heigth of guide exit (h2) */ div1 = ((lbh/2-louth)*(lbh/2-louth))/a_ell_q; h2 = sqrt(b_ell_q*(1-div1)); h2 = h2*2; } } else if (!strcmp(fu,"parabolical")) { if ((linh > 0) && (louth > 0)) { fprintf(stderr,"Component: %s (Guide_tapering) Two focal\n",NAME_CURRENT_COMP); fprintf(stderr," points lout and linh are not allowed! \n", file_name); free(fu);exit(-1); } if (louth == 0 && linh == 0) { /* plane mirrors (vertical) */ h2 = h1; } else { /* parabolical mirrors */ if (linh == 0) { Div1=((2*louth+2*l)*(2*louth+2*l))/4; p2_para=((sqrt(Div1+(h12*h12)))-(louth+l))*2; /* calculate heigth of guide exit (h2) */ h2 = sqrt(p2_para*(louth+p2_para/4)); h2 = h2*2; } else { /* anti-trompete */ Div1=((2*linh)*(2*linh))/4; p2_para=((sqrt(Div1+(h12*h12)))-linh)*2; /* calculate heigth of guide exit (h2) */ h2 = sqrt(p2_para*(l+linh+p2_para/4)); h2 = h2*2; } } } else if (!strncmp(fu,"file",4)) { pos = strtok(fu,"="); while (pos=strtok(0,"=")) { strcpy(file_name,pos); } if ((num=fopen(file_name,"r")) == NULL) { fprintf(stderr,"Component: %s (Guide_tapering)\n",NAME_CURRENT_COMP); fprintf(stderr," File %s not found! \n", file_name); free(fu);exit(-1); } else { ii = 0; while (!feof(num)) { fgets(segment[ii].st,128,num); if (ii > 799) { fprintf(stderr,"Number of segments is limited to 800 !! \n",NAME_CURRENT_COMP); free(fu);exit(-1); } ii++; } fclose(num); ii--; } seg = ii-3; l_seg=l/seg; for (i=3;i 0)) { for (i=1;i<(seg+1);i++) { h1_in[i] = (sqrt((p2_para/4+louth+(l_seg*ii))*p2_para))*2; ii=ii-1; h2_out[i-1] = h1_in[i]; } } else { for (i=1;i<(seg+1);i++) { h1_in[i] = (sqrt((p2_para/4+linh+(l_seg*i))*p2_para))*2; h2_out[i-1] = h1_in[i]; } } } } /* compute each value for horizontal mirrors */ w12 = w1/2; if (!strcmp(fu,"elliptical")) { /* calculate lbw the distance between focal points of horizontal mirrors */ lbw = l + linw + loutw; /* calculate parameter b of elliptical equestion - horizontal mirrors */ if (linw == 0 && loutw == 0 ) { /* plane mirrors (horizontal) */ b_ell_q = 0; w2 = w1; } else { /* elliptical mirrors */ u1 = ((lbw/2)-linw)*((lbw/2)-linw); u2 = (lbw/2)*(lbw/2); div1 = u1/u2; b_ell_q = (w12*w12)/(1-div1); u1 = sqrt((linw*linw)+(w12*w12)); u2 = sqrt((w12*w12) + ((l+loutw)*(l+loutw))); a_ell_q = ((u1 + u2)/2)*((u1 + u2)/2); b_ell_q = a_ell_q - ((lbw/2)*(lbw/2)); /* calculate weigth of guide exit (w2) */ div1 = ((lbw/2-loutw)*(lbw/2-loutw))/a_ell_q; w2 = sqrt(b_ell_q*(1-div1)); w2 = w2*2; } } else if (!strcmp(fu,"parabolical")) { if ((linw > 0) && (loutw > 0)) { fprintf(stderr,"Component: %s (Guide_tapering) Two focal\n",NAME_CURRENT_COMP); fprintf(stderr," points linw and loutw are not allowed! \n", file_name); free(fu);exit(-1); } if (loutw == 0 && linw == 0) { /* plane mirrors (horizontal) */ w2 = w1; } else { if (linw == 0) { /* parabolical mirrors */ Div1=((2*loutw+2*l)*(2*loutw+2*l))/4; p2_para=((sqrt(Div1+(w12*w12)))-(loutw+l))*2; /* calculate weigth of guide exit (w2) */ w2 = sqrt(p2_para*(loutw+p2_para/4)); w2 = w2*2; } else { /* anti-trompete */ Div1=((2*linw)*(2*linw))/4; p2_para=((sqrt(Div1+(w12*w12)))-linw)*2; /* calculate heigth of guide exit (w2) */ w2 = sqrt(p2_para*(l+linw+p2_para/4)); w2 = w2*2; } } } fprintf(stderr,"Component: %s (Guide_tapering)\n",NAME_CURRENT_COMP); fprintf(stderr," Weigth at the guide exit (w2): %lf \n", w2); if (w2 <= 0) { fprintf(stderr,"Component: %s (Guide_tapering)\n", NAME_CURRENT_COMP); fprintf(stderr," Weigth at the guide exit (w2) was calculated\n"); fprintf(stderr," <=0; Please change the parameter w1 and/or\n"); fprintf(stderr," l! \n"); free(fu);exit(-1); } if (!strcmp(fu,"elliptical")) { w1_in[0]=w1; a_ell_q = (lbw/2)*(lbw/2); for (i=1;i 0)) { for (i=1;i<(seg+1);i++) { w1_in[i] = (sqrt((p2_para/4+loutw+(l_seg*ii))*p2_para))*2; ii=ii-1; w2_out[i-1] = w1_in[i]; } } else { for (i=1;i<(seg+1);i++) { w1_in[i] = (sqrt((p2_para/4+linw+(l_seg*i))*p2_para))*2; w2_out[i-1] = w1_in[i]; } } } } free(fu); for (i=0;i= w1_in[ii]/2.0 || y <= -hhalf[ii] || y >= hhalf[ii]) ABSORB; /* Shift origin to center of channel hit (absorb if hit dividing walls) */ x += w1_in[ii]/2.0; edge = floor(x/w1c[ii])*w1c[ii]; if(x - edge > w1c[ii]) { x -= w1_in[ii]/2.0; /* Re-adjust origin */ ABSORB; } x -= (edge + (w1c[ii]/2.0)); hadj = edge + (w1c[ii]/2.0) - w1_in[ii]/2.0; for(;;) { /* Compute the dot products of v and n for the four mirrors. */ av = l_seg*vx; bv = ww[ii]*vz; ah = l_seg*vy; bh = hh[ii]*vz; vdotn_v1 = bv + av; /* Left vertical */ vdotn_v2 = bv - av; /* Right vertical */ vdotn_h1 = bh + ah; /* Lower horizontal */ vdotn_h2 = bh - ah; /* Upper horizontal */ /* Compute the dot products of (O - r) and n as c1+c2 and c1-c2 */ cv1 = -whalf[ii]*l_seg - z*ww[ii]; cv2 = x*l_seg; ch1 = -hhalf[ii]*l_seg - z*hh[ii]; ch2 = y*l_seg; /* Compute intersection times. */ t1 = (l_seg - z)/vz; i = 0; if(vdotn_v1 < 0 && (t2 = (cv1 - cv2)/vdotn_v1) < t1) { t1 = t2; i = 1; } if(vdotn_v2 < 0 && (t2 = (cv1 + cv2)/vdotn_v2) < t1) { t1 = t2; i = 2; } if(vdotn_h1 < 0 && (t2 = (ch1 - ch2)/vdotn_h1) < t1) { t1 = t2; i = 3; } if(vdotn_h2 < 0 && (t2 = (ch1 + ch2)/vdotn_h2) < t1) { t1 = t2; i = 4; } if(i == 0) { break; /* Neutron left guide. */ } PROP_DT(t1); switch(i) { case 1: /* Left vertical mirror */ nlen2 = l_seg*l_seg + ww[ii]*ww[ii]; q = V2Q*(-2)*vdotn_v1/sqrt(nlen2); dd = 2*vdotn_v1/nlen2; vx = vx - dd*l_seg; vz = vz - dd*ww[ii]; break; case 2: /* Right vertical mirror */ nlen2 = l_seg*l_seg + ww[ii]*ww[ii]; q = V2Q*(-2)*vdotn_v2/sqrt(nlen2); dd = 2*vdotn_v2/nlen2; vx = vx + dd*l_seg; vz = vz - dd*ww[ii]; break; case 3: /* Lower horizontal mirror */ nlen2 = l_seg*l_seg + hh[ii]*hh[ii]; q = V2Q*(-2)*vdotn_h1/sqrt(nlen2); dd = 2*vdotn_h1/nlen2; vy = vy - dd*l_seg; vz = vz - dd*hh[ii]; break; case 4: /* Upper horizontal mirror */ nlen2 = l_seg*l_seg + hh[ii]*hh[ii]; q = V2Q*(-2)*vdotn_h2/sqrt(nlen2); dd = 2*vdotn_h2/nlen2; vy = vy + dd*l_seg; vz = vz - dd*hh[ii]; break; } /* Now compute reflectivity. */ if((i <= 2 && mx == 0) || (i > 2 && my == 0)) { x += hadj; /* Re-adjust origin */ ABSORB; } if((i <= 2 && q > Qcx) || (i > 2 && q > Qcy)) { if (i <= 2) { double arg = (q - mx*Qcx)/W; if(arg < 10) p *= .5*(1-tanh(arg))*(1-alphax*(q-Qcx)); else { x += hadj; /* Re-adjust origin */ ABSORB; /* Cutoff ~ 1E-10 */ } } else { double arg = (q - my*Qcy)/W; if(arg < 10) p *= .5*(1-tanh(arg))*(1-alphay*(q-Qcy)); else { x += hadj; /* Re-adjust origin */ ABSORB; /* Cutoff ~ 1E-10 */ } } } p *= R0; x += hadj; SCATTER; x -= hadj; } x += hadj; /* Re-adjust origin */ /* store current neutron state */ if (ii<(seg-1)) { z=z-l_seg; STORE_NEUTRON(INDEX_CURRENT_COMP,x,y,z,vx,vy,vz,t,sx,sy,sz,p); } } if (seg > 1) { z=z+l_seg*(seg-1); } %} MCDISPLAY %{ double x; int i,ii; magnify("xy"); for (ii=0; ii < segno; ii++) { for(i = 0; i < 1; i++) { multiline(5, i*w1c[ii] - w1_in[ii]/2.0, -h1_in[ii]/2.0, (double)l_seg*ii, i*w2c[ii] - w2_out[ii]/2.0, -h2_out[ii]/2.0, (double)l_seg*(ii+1.0), i*w2c[ii] - w2_out[ii]/2.0, h2_out[ii]/2.0, (double)l_seg*(ii+1.0), i*w1c[ii] - w1_in[ii]/2.0, h1_in[ii]/2.0, (double)l_seg*ii, i*w1c[ii] - w1_in[ii]/2.0, -h1_in[ii]/2.0, (double)l_seg*ii); multiline(5, (i+1.0)*w1c[ii] - w1_in[ii]/2.0, -h1_in[ii]/2.0, (double)l_seg*ii, (i+1.0)*w2c[ii] - w2_out[ii]/2.0, -h2_out[ii]/2.0, (double)l_seg*(ii+1.0), (i+1.0)*w2c[ii] - w2_out[ii]/2.0, h2_out[ii]/2.0, (double)l_seg*(ii+1.0), (i+1.0)*w1c[ii] - w1_in[ii]/2.0, h1_in[ii]/2.0, (double)l_seg*ii, (i+1.0)*w1c[ii] - w1_in[ii]/2.0, -h1_in[ii]/2.0, (double)l_seg*ii); } } line(-w1/2.0, -h1/2.0, 0.0, w1/2.0, -h1/2.0, 0.0); for(i=1; i