/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright (C) 1997-2008, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: StatisticalChopper * * %I * Written by: C. Monzat/E. Farhi/S. Rozenkranz * Date: Nov 10th, 2009 * Origin: ILL * * Statistical (correlation) Chopper * * %D * This component is a statistical chopper based on a pseudo random sequence that the user * can specify. Inspired from DiskChopper, it should be used with SatisticalChopper_Monitor component. * * Example: StatisticalChopper(nu=1487*60/255) use default sequence * * %P * INPUT PARAMETERS: * * yheight: [m] Slit height (if = 0, equal to radius). Auto centering of beam at half height. * radius: [m] Radius of the disc * nu: [Hz] Frequency of the Chopper, omega=2*PI*nu * (algebraic sign defines the direction of rotation) * sequence: [str] Slit sequence around disk, with 0=closed, 1=open positions. * NULL or "" will set a default sequence. * * Optional parameters: * isfirst: [0/1] Set it to 1 for the first chopper position in a cw source * (it then spreads the neutron time distribution) * n_pulse: [1] Number of pulses (Only if isfirst) * jitter: [s] Jitter in the time phase * abs_out: [0/1] Absorb neutrons hitting outside of chopper radius? * delay: [s] Time 'delay'. * phase: [deg] Angular 'delay' (overrides delay) * verbose: [1] Set to 1 to display Disk chopper configuration * * %L * R. Von Jan and R. Scherm. The statistical chopper for neutron time-of-flight spectroscopy. Nuclear Instruments and Methods, 80 (1970) 69-76. * * %E *******************************************************************************/ DEFINE COMPONENT StatisticalChopper SETTING PARAMETERS (radius=0.5, yheight=0, nu, jitter=0, delay=0, isfirst=0, abs_out=1, phase=0, verbose=0, n_pulse=1, string sequence="NULL") OUTPUT PARAMETERS () /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ DECLARE %{ double delta_y; double height; double omega; int *Sequence; /* Array containing 0=closed and 1=opened*/ int *SequenceIndex; /* index of '1' in the sequence */ int nslit; int m; /* sum of opened slits ('1' in sequence) */ %} INITIALIZE %{ Sequence=NULL; SequenceIndex=NULL; nslit=0; m=0; /* If slit height 'unset', assume full opening */ if (yheight == 0) { height=radius; } else { height=yheight; } delta_y = radius-height/2; /* radius at beam center */ omega=2.0*PI*nu; /* rad/s */ if (!sequence) { fprintf(stderr,"StatisticalChopper: %s: sequence is NULL. \n", NAME_CURRENT_COMP); exit(-1); } if (strlen(sequence) <=2) { fprintf(stderr,"StatisticalChopper: %s: Sequence is too short (length=%i). Use conventional DiskChopper instead.\n", NAME_CURRENT_COMP, nslit); exit(-1); } if (!strlen(sequence) || !strcmp(sequence,"NULL")) strcpy(sequence, "100000010101001101110010010011110111111001011000011101001010111100101011111011010011011111000100011011100010000010001101101010100001110111001011111011001110001100011011011110101010000000110000011110011010100101011100000001001101000111100010110110001100100" ); if (yheight && yheight>radius) { fprintf(stderr,"StatisticalChopper: %s: yheight must be < radius\n", NAME_CURRENT_COMP); exit(-1); } if (isfirst && n_pulse <=0) { fprintf(stderr,"StatisticalChopper: %s: wrong First chopper pulse number (n_pulse=%g)\n", NAME_CURRENT_COMP, n_pulse); exit(-1); } if (!omega) { fprintf(stderr,"StatisticalChopper: %s WARNING: chopper frequency is 0!\n", NAME_CURRENT_COMP); omega = 1e-15; /* We should actually use machine epsilon here... */ } if (!abs_out) { fprintf(stderr,"StatisticalChopper: %s WARNING: chopper will NOT absorb neutrons outside radius %g [m]\n", NAME_CURRENT_COMP, radius); } /* Calulate delay from phase and vice versa, 'direction' moderated by sign of omega */ delay *=omega/fabs(omega); if (phase) { if (delay) { fprintf(stderr,"StatisticalChopper: %s WARNING: delay AND phase specified. Using phase setting\n", NAME_CURRENT_COMP); } phase*=DEG2RAD; /* 'Delay' should always be a delay, taking rotation direction into account: */ delay=omega*phase/(omega*omega); } else { phase=delay*omega; /* rad */ } if (verbose && nu) { printf("StatisticalChopper: %s: frequency=%g [Hz] %g [rpm] time frame=%g [s] phase=%g [deg]\n", NAME_CURRENT_COMP, nu, nu*60, fabs(1/nu), phase*RAD2DEG); printf(" height=%g [m], slits centered at radius=%g [m]\n", height, delta_y); } nslit=strlen(sequence); if (nslit>1) { int i=0,j=0; double c=0; int index=0; Sequence = malloc(nslit * sizeof(int)); SequenceIndex = malloc(nslit * sizeof(int)); if (!Sequence){ fprintf(stderr,"StatisticalChopper: %s: Memory exhausted when allocating chopper sequence.\n", NAME_CURRENT_COMP); exit(-1); } if (!SequenceIndex){ fprintf(stderr,"StatisticalChopper: %s: Memory exhausted when allocating chopper sequence table lookup.\n", NAME_CURRENT_COMP); } /* build sequence index (where sequence==1) */ for (i=0;iradius*radius) { ABSORB; } /* Does neutron hit inner solid part of chopper in case of yheight!=radius ? */ if ((x*x+yprime*yprime)<(radius-height)*(radius-height)) { ABSORB; } if (isfirst && SequenceIndex) { /* choose a slit in the sequence and set time accordingly */ t = atan2(x,yprime)/omega + SequenceIndex[(int)(rand01()*m)]/nslit/nu -delay +jitter*randnorm()+ (n_pulse > 1 ? floor(n_pulse*rand01()/fabs(nu)) : 0); p *= (double)m/(double)nslit; /* transmission */ } else { /* check if we pass through a slit */ int seq; double angle = 2*PI*nu*(t-atan2(x,yprime)/omega+jitter*randnorm()) - phase; seq=(int) floor(angle*nslit/(2*PI)) % nslit; if ( seq < 0 ) seq += nslit; if (Sequence[seq]) SCATTER; else ABSORB; } %} FINALLY %{ free(Sequence); Sequence=NULL; free(SequenceIndex); SequenceIndex=NULL; %} MCDISPLAY %{ int i; double radius_min=radius-height; circle("xy", 0.0, -radius, 0.0, radius); for (i=0; i