From 9865c3f9c5e431227927f1966545a76c05d3cb41 Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Fri, 11 Oct 2024 12:34:12 +0200
Subject: [PATCH] Initial footwork toward 
 https://github.com/McStasMcXtrace/McCode/issues/1736

---
 mcxtrace-comps/optics/Collimator_linear.comp | 120 +++++++++++++++++++
 1 file changed, 120 insertions(+)
 create mode 100644 mcxtrace-comps/optics/Collimator_linear.comp

diff --git a/mcxtrace-comps/optics/Collimator_linear.comp b/mcxtrace-comps/optics/Collimator_linear.comp
new file mode 100644
index 000000000..aad860329
--- /dev/null
+++ b/mcxtrace-comps/optics/Collimator_linear.comp
@@ -0,0 +1,120 @@
+/*******************************************************************************
+*
+* McXtrace, x-ray tracing package
+*         Copyright, All rights reserved
+*         DTU Physics, Kgs. Lyngby, Denmark
+*         Synchrotron SOLEIL, Saint-Aubin, France
+*
+* Component: Collimator_linear
+*
+* %I
+* Written by: Peter Willendrup
+* Date: October 2024
+* Origin: DTU
+* Modified by: Peter Willendrup from McStas Collimator_linear from Krisitan Nielsen
+*
+* A simple analytical Soller collimator (with triangular transmission).
+*
+* %D
+* Soller collimator with rectangular opening and specified length. The
+* transmission function is an average and does not utilize knowledge of the
+* actual neutron trajectory. A zero divergence disables collimation (then the
+* component works as a double slit).
+*
+* Example: Collimator_linear(xmin=-0.1, xmax=0.1, ymin=-0.1, ymax=0.1, length=0.25, divergence=40,transmission=0.7)
+*
+* %P
+* INPUT PARAMETERS:
+*
+* xmin: [m]                      Lower x bound on slits
+* xmax: [m]                      Upper x bound on slits
+* ymin: [m]                      Lower y bound on slits
+* ymax: [m]                      Upper y bound on slits
+* xwidth: [m]                    Width of slits
+* yheight: [m]                   Height of slits
+* length: [m]                    Distance between input and output slits
+* divergence: [minutes of arc]   Divergence horizontal angle (calculated as atan(d/length), where d is the blade spacing)
+* transmission: [1]              Transmission of Soller (0<=t<=1)
+* divergenceV: [minutes of arc]  Divergence vertical angle
+*
+* %E
+*******************************************************************************/
+DEFINE COMPONENT Collimator_linear
+DEFINITION PARAMETERS ()
+SETTING PARAMETERS (xmin=-0.01, xmax=0.01, ymin=-0.025, ymax=0.025, xwidth=0, yheight=0, length=0.1, divergence=10, transmission=1, divergenceV=0)
+OUTPUT PARAMETERS ()
+/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */
+
+DECLARE
+%{
+  double slope;
+  double slopeV;
+%}
+INITIALIZE
+%{
+  slope = tan(MIN2RAD*divergence);
+  slopeV= tan(MIN2RAD*divergenceV);
+  if (xwidth  > 0) { xmax = xwidth/2;  xmin = -xmax; }
+  if (yheight > 0) { ymax = yheight/2; ymin = -ymax; }
+
+  if ((xmin >= xmax) || (ymin >= ymax)) {
+    printf("Collimator_linear: %s: Null slit opening area !\n"
+	         "ERROR              (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting",
+           NAME_CURRENT_COMP);
+    exit(0);
+  }
+%}
+TRACE
+%{
+    double Phi, dl;
+
+    PROP_Z0;
+    if (x<xmin || x>xmax || y<ymin || y>ymax)
+      ABSORB;
+
+    double dirx,diry,dirz;
+
+    /* Use direction vector z-component to decide needed propagation length */
+    dirx=kx;diry=ky;dirz=kz;
+    NORM(dirx,diry,dirz);
+   
+    PROP_DL(length/dirz);
+    if (x<xmin || x>xmax || y<ymin || y>ymax)
+      ABSORB;
+
+    if(slope > 0.0)
+    {
+      Phi = fabs(kx/kz);
+      if (Phi > slope)
+        ABSORB;
+      else
+        p *= transmission*(1.0 - Phi/slope);
+      SCATTER;
+    }
+    if (slopeV > 0) {
+      Phi = fabs(ky/kz);
+      if (Phi > slopeV)
+        ABSORB;
+      else
+        p *= transmission*(1.0 - Phi/slopeV);
+      SCATTER;
+    }
+%}
+
+MCDISPLAY
+%{
+  double x;
+  int i;
+
+  
+  for(x = xmin, i = 0; i <= 10; i++, x += (xmax - xmin)/10.0)
+    multiline(5, x, (double)ymin, 0.0, x, (double)ymax, 0.0,
+              x, (double)ymax, (double)length, x, (double)ymin, (double)length,
+              x, (double)ymin, 0.0);
+  line(xmin, ymin, 0,   xmax, ymin, 0);
+  line(xmin, ymax, 0,   xmax, ymax, 0);
+  line(xmin, ymin, length, xmax, ymin, length);
+  line(xmin, ymax, length, xmax, ymax, length);
+%}
+
+END