Initial commit of Phonon_BvK_PG and test instrument

mcstas-comps/contrib/Phonon_BvK_PG.comp

@@ -38,7 +38,7 @@
 * yheight: [m]        Height of sample in y direction
 * xwidth: [m]         Horiz. dimension of sample if slap
 * zdepth: [m]         Depth of sample if slap
-
+*
 * sigma_abs: [barns]  Absorption cross section at 2200 m/s per atom
 * sigma_inc: [barns]  Incoherent scattering cross section per atom
 * DW: [1]             Debye-Waller factor; scalar value (no q-dependence)
@@ -71,7 +71,6 @@ DEFINE COMPONENT Phonon_BvK_PG
 SETTING PARAMETERS (hh=0,kk=0,ll=0,radius=0, xwidth=0, yheight=0, zdepth=0, sigma_abs=0,sigma_inc=0,DW=1,T,
 focus_r=0,focus_xw=0,focus_yh=0,focus_aw=0,focus_ah=0,target_x=0, target_y=0, target_z=0, int target_index=0, int mode_input, int e_steps_low, int e_steps_high, int verbose_input=0, int dispersion=0)
 NOACC
-// The component is currently "NOACC" only.
 
 SHARE
 %{
@@ -1121,12 +1120,6 @@ MCDISPLAY
   } else if (xwidth && yheight) { /* box/rectangle */
     box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0);
   }
-  //  circle("xz", 0,  yheight/2.0, 0, radius);
-  //  circle("xz", 0, -yheight/2.0, 0, radius);
-  //  line(-radius, -yheight/2.0, 0, -radius, +yheight/2.0, 0);
-  //  line(+radius, -yheight/2.0, 0, +radius, +yheight/2.0, 0);
-  //  line(0, -yheight/2.0, -radius, 0, +yheight/2.0, -radius);
-  //  line(0, -yheight/2.0, +radius, 0, +yheight/2.0, +radius);
 %}
 
 END

mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/Test_Phonon_BvK_PG.instr

@@ -0,0 +1,315 @@
+/*******************************************************************************
+* Instrument: Test_Phonon_BvK
+*
+* %I
+* Written by: Kim Lefmann
+* Date: 20.07.2025
+* Origin: NBI, KU
+* %INSTRUMENT_SITE: Generic
+*
+* Simple thermal triple axis spectrometer for test of Phonon_BvK_PG component
+* Modified from older version by Dorothy Wang, Kim Lefmann
+*
+* %D
+* Simple TAS to test the Phonon_BvK_PG component
+*
+* %Example: -y Detector: e_mon_beforeana_zoom_I=6.90143e-16
+*
+* %P
+* E:       [meV] Energy transfer
+* Ef:      [meV] Final energy
+* Dlambda: [AA] width of wavelength band
+* h:       [rlu] q-value in hexagonal plane
+* l:       [rlu] q-value out of hexagonal plane
+* dA3:     [deg] A3 offset
+* phononmode: [1] choice of a specific phonon mode  
+*
+* %L
+* <reference/HTML link>
+*
+* %E
+*******************************************************************************/
+DEFINE INSTRUMENT Test_Phonon_BvK_PG(E=4.5,Ef=25,Dlambda=0.05,h=0,l=3.5,dA3=90, Temp=2, width = 0.001, coll = 20, int phononmode=0,int E_steps_high=100, int E_steps_low=100, int Verbose=0, int DISP=0)
+
+//SHELL " cc -c eigen.c "
+//DEPENDENCY " eigen.o "
+
+DECLARE
+%{
+  double Gqx,Gqy,Gqz;
+//double Ef=24.8
+double Ei;
+//meV
+double qx,qy,qz;
+double thetaM;
+double twothetaS;
+double thetaA;
+double A3;
+//double pi;
+double QM;
+double alpha;
+double lambda_i;
+double SMALL__NUMBER;
+%}
+
+INITIALIZE
+%{
+// Set monochromator/analyzer Q-value for PG
+QM = 1.8734;
+SMALL__NUMBER = 1e-6;
+//pi = 3.14159265;
+
+double a = 2.461;
+double c =  6.708;
+//PG lattice constants, in AA, same as in Phonon_BvK_PG.comp
+
+double astar = 4*pi/(sqrt(3)*a);
+double cstar = 2*pi/c;
+//PG reciprocal lattice vectors, in 1/AA
+
+//calculate Ei
+Ei=Ef+E;
+
+//calculate ki, kf, lambda_i, q
+double ki = V2K*SE2V*sqrt(Ei);
+double kf = V2K*SE2V*sqrt(Ef);
+lambda_i=2*pi/ki;
+double q = sqrt(h*h*astar*astar+l*l*cstar*cstar);
+
+//calculate 2thetaM and 2thetaA
+thetaM = asin(QM/(2*ki))*RAD2DEG;
+thetaA = asin(QM/(2*kf))*RAD2DEG;
+
+//calculate scattering angle and sample rotation
+twothetaS = acos((q*q-ki*ki-kf*kf)/(-2*ki*kf))*RAD2DEG;
+alpha = acos((kf*kf-ki*ki-q*q)/(-2*ki*q));
+A3=(asin(l*cstar/(q+SMALL__NUMBER))-alpha)*RAD2DEG;
+
+//Printout calculations
+printf("a = %g c = %g astar = %g cstar = %g \n",a,c,astar,cstar);
+printf("ki = %g kf = %g q = %g lambda_i = %g\n",ki,kf,q,lambda_i);
+printf("thetaA = %g thetaM = %g\n",thetaA,thetaM);
+printf("twothetaS = %g \n",twothetaS);
+printf("A3 = %g \n",A3);
+printf("alpha = %g \n",alpha);
+
+%}
+
+TRACE
+
+COMPONENT origin = Progress_bar()
+AT (0, 0, 0) RELATIVE ABSOLUTE
+
+COMPONENT source = Source_Maxwell_3(
+    xwidth=width,
+    yheight=width,
+    Lmin=lambda_i-Dlambda/2, 
+    Lmax=lambda_i+Dlambda/2, 
+    dist=7.5, 
+    focus_yh=width, 
+    focus_xw=width, 
+    T1=300,     T2=300,     T3=300, 
+    I1=1E15,    I2=1E15,    I3=1E15)
+AT (0, 0, 0) RELATIVE PREVIOUS
+
+COMPONENT l_monitor_before_mono = L_monitor(
+    nL=200, 
+    filename="L_before_mono", 
+    Lmin=0, 
+    Lmax=4,   
+    xwidth=0.16, 
+    yheight=0.25, 
+    restore_neutron=1)
+AT (0, 0, 7) RELATIVE PREVIOUS
+
+COMPONENT monochromator_flat = Monochromator_flat(
+     mosaicv=30,
+     mosaich=30,
+     zwidth = width,
+     yheight = width,
+     Q=QM)
+AT (0, 0, 7.5) RELATIVE source
+ROTATED (0, thetaM, 0) RELATIVE source
+
+COMPONENT arm1 = Arm()
+AT (0, 0, 0) RELATIVE PREVIOUS
+ROTATED (0, thetaM, 0) RELATIVE PREVIOUS
+
+COMPONENT collimator_linear1 = Collimator_linear(
+    xwidth=0.1, 
+    yheight=0.2, 
+    length=0.2, 
+    divergence=coll)
+AT (0, 0, 0.6) RELATIVE arm1
+
+COMPONENT l_monitor_before_sample_big = L_monitor(
+    nL=300, 
+    filename="L_before_sample", 
+    Lmin=lambda_i - Dlambda, 
+    Lmax=lambda_i + Dlambda,
+    xwidth = 0.1,
+    yheight = 0.1,
+    restore_neutron=1)
+AT (0, 0, 0.9) RELATIVE arm1
+
+COMPONENT e_monitor_before_sample_big = E_monitor(
+    nE=100, 
+    filename="E_before_sample", 
+    Emin=Ei-3, 
+    Emax=Ei+3, 
+    xwidth = 0.1,
+    yheight = 0.1,
+   restore_neutron=1)
+AT (0, 0, 0.01) RELATIVE PREVIOUS
+
+COMPONENT PSD_monitor_before_sample = PSD_monitor(
+    nx = 200,
+    ny = 200,
+    xwidth = width/2,
+    yheight = width/2,
+    filename="PSD_before_sample", 
+    restore_neutron=1)
+AT (0, 0, 0.998) RELATIVE arm1
+
+COMPONENT arm2 = Arm()
+AT (0, 0, 1) RELATIVE arm1
+ROTATED (0, -A3+dA3, 0) RELATIVE arm1
+
+SPLIT 10
+/*COMPONENT incoherent = Incoherent( 
+    yheight=width/2, 
+    xwidth = width/2,
+   zdepth = width/2000,
+    focus_xw=width, 
+    focus_yh=width, 
+    target_index=3,
+   order = 1, 
+    sigma_abs=0) */
+COMPONENT phonon_bvk_pg = Phonon_BvK_PG(
+//COMPONENT phonon_bvk_pg = Phonon_BvK_PG_eigenvector_Dec2024(
+//    radius=width/2,
+    xwidth=width/2, zdepth=width/20,
+    yheight=width/2, 
+    sigma_abs=0, 
+    sigma_inc=0, 
+    DW=1, 
+    T=Temp, 
+    mode_input=phononmode,
+    target_index=3,
+    focus_xw=width,
+    focus_yh=width,
+    e_steps_low = E_steps_low,
+    e_steps_high = E_steps_high,
+    verbose_input=Verbose,
+    hh = h, kk = 0, ll = l, dispersion = DISP
+)
+AT (0, 0, 0) RELATIVE PREVIOUS
+ROTATED (0, 0, 0) RELATIVE PREVIOUS  // To use the (0,k,l) scattering plane, rotate (0, 0, -60); for (h,k,0) rotate (90, 0, 0)
+
+COMPONENT arm3 = Arm()
+AT (0, 0, 0) RELATIVE PREVIOUS
+ROTATED (0, -twothetaS, 0) RELATIVE arm1
+
+/*COMPONENT PSD_sample_image = PSD_monitor(
+    nx=200, 
+    ny=200, 
+   xwidth = 0.01,
+   yheight = 0.01,
+    filename="PSD_sample_image", 
+    restore_neutron=1)
+AT (0, 0, 1.005) RELATIVE arm1*/
+
+/*COMPONENT collimator_linear2 = Collimator_linear(
+    xwidth=0.1, 
+    yheight=0.2, 
+    length=0.2, 
+    divergence=coll)
+AT (0, 0, 0.5) RELATIVE PREVIOUS
+*/
+
+COMPONENT PSD_monitor_after_sample = PSD_monitor(
+    nx=200, 
+    ny=200, 
+    filename="PSD_after_sample", 
+    restore_neutron=1)
+AT (0, 0, 0.9999) RELATIVE arm3
+
+COMPONENT slit1 = Slit(
+    xwidth=width, 
+    yheight=width)
+AT (0, 0, 1) RELATIVE arm3
+
+COMPONENT e_monitor_beforeana = E_monitor(
+    nE=750,
+    filename="E_before_ana", 
+    xwidth=2*width, 
+    yheight=2*width, 
+    Emin=0, 
+    Emax=75)
+AT (0, 0, 0.001) RELATIVE PREVIOUS
+
+COMPONENT e_mon_beforeana_zoom = E_monitor(
+    nE=750,
+    filename="E_before_ana_zoom", 
+    xwidth=2*width, 
+    yheight=2*width, 
+    Emin=Ef*0.9, 
+    Emax=Ef*1.1)
+AT (0, 0, 0.001) RELATIVE PREVIOUS
+COMPONENT analyzer = Monochromator_flat(
+    zwidth=4*width,
+    yheight=4*width,
+    mosaicv=30,
+    mosaich=30,
+    Q=QM
+    )
+AT (0, 0, 0.1) RELATIVE PREVIOUS
+ROTATED (0, thetaA, 0) RELATIVE PREVIOUS
+
+COMPONENT arm4 = Arm()
+AT (0, 0, 0) RELATIVE PREVIOUS
+ROTATED (0, thetaA, 0) RELATIVE PREVIOUS
+
+COMPONENT slit2 = Slit(
+    xwidth=width,
+    yheight=width)
+AT (0, 0, 1) RELATIVE PREVIOUS
+EXTEND
+%{
+//    qx = MC_GETPAR(phonon_bvk_pg,q_x);
+//    qy = MC_GETPAR(phonon_bvk_bg,q_y);
+//    qz = MC_GETPAR(phonon_bvk_bg,q_z);
+//    printf("hit; q = (%g, %g, %g)\n",0,0,0);
+%}
+
+COMPONENT Emon_after_analyzer = E_monitor(
+    nE=750,
+    filename="E_after_analyzer", 
+    xwidth=0.05,
+    yheight=0.10, 
+    Emin=0, 
+    Emax=250)
+AT (0, 0, 0.001) RELATIVE PREVIOUS
+
+COMPONENT Emon_1storder = E_monitor(
+    nE=50,
+    filename="E_1st_order", 
+    xwidth=width,
+    yheight=width, 
+    Emin=Ef*0.9, 
+    Emax=Ef*1.1)
+AT (0, 0, 0.001) RELATIVE PREVIOUS
+
+/*COMPONENT PSD_detector = PSD_monitor(
+    nx=128, 
+    ny=128, 
+    filename="PSD_detector", 
+    xwidth=0.05,
+    yheight=0.10)
+AT (0, 0, 0.001) RELATIVE PREVIOUS
+*/
+FINALLY
+%{
+%}
+
+END