From 6b84b262513a94e26855f2a927511b278df9bdb8 Mon Sep 17 00:00:00 2001 From: Kim Lefmann Date: Tue, 24 Feb 2026 11:56:33 +0100 Subject: [PATCH 01/15] Initial commit of Phonon_BvK_PG and test instrument --- mcstas-comps/contrib/Phonon_BvK_PG.comp | 1334 +++++++++++++++++ .../Test_Phonon_BvK_PG.instr | 315 ++++ .../Tests_samples/Test_Phonon_BvK_PG/eigen.c | 451 ++++++ .../Tests_samples/Test_Phonon_BvK_PG/eigen.h | 40 + .../Test_Phonon_BvK_PG/mergesort.c | 73 + .../Tests_samples/Test_Phonon_BvK_PG/types.h | 50 + 6 files changed, 2263 insertions(+) create mode 100644 mcstas-comps/contrib/Phonon_BvK_PG.comp create mode 100644 mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/Test_Phonon_BvK_PG.instr create mode 100644 mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/eigen.c create mode 100644 mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/eigen.h create mode 100644 mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/mergesort.c create mode 100644 mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/types.h diff --git a/mcstas-comps/contrib/Phonon_BvK_PG.comp b/mcstas-comps/contrib/Phonon_BvK_PG.comp new file mode 100644 index 000000000..c2e93cdb4 --- /dev/null +++ b/mcstas-comps/contrib/Phonon_BvK_PG.comp @@ -0,0 +1,1334 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2000 Risoe National Laboratory. +* +* %I +* Written by: Kim Lefmann +* Date: 06.12.2025 +* Origin: NBI, KU +* +* A sample for phonon scattering based on cross section expressions from the Boothroyd textbook +* based upon the algorithm from component Phonon_simple (with expressions from Squires, Ch. 3) +* Using the Born-von Karman force constnts and the normal mode formalism +* +* Optimized for performance from the previous Phonon_BvK_PG_eígenvector:Dec2024.comp and Phonon_BvK_PG_eígenvector:July2025.comp +* +* %D +* Single-cylinder or slap shaped shape. +* Absorption included. +* No multiple scattering. +* No incoherent scattering emitted, but incoherent is present in attenuation +* No attenuation from coherent scattering. No Bragg scattering. +* Specialized for PG +* +* Algorithm: +* 0. Always perform the scattering if possible (otherwise ABSORB) +* 1. Choose direction within a focusing solid angle +* 2. Select a phonon mode at random (alternatively mode number fixed by user) +* 3. Calculate the zeros of (E_i-E_f-hbar omega(kappa)) as a function of k_f +* 4. Choose one value of k_f (there always at least one possible; see e.g. Squires) +* 5. Perform the correct weight transformation +* For details: see article A. F. Davidsen et al, manuscript for J Neutron Res 2025 +* +* %P +* INPUT PARAMETERS: +* radius: [m] Outer radius of sample in (x,z) plane +* 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) +* T: [K] Temperature +* focus_r: [m] Radius of sphere containing target. +* focus_xw: [m] horiz. dimension of a rectangular area +* focus_yh: [m] vert. dimension of a rectangular area +* focus_aw: [deg] horiz. angular dimension of a rectangular area +* focus_ah: [deg] vert. angular dimension of a rectangular area +* target_x: [m] position of target to focus at . Transverse coordinate +* target_y: [m] position of target to focus at. Vertical coordinate +* target_z: [m] position of target to focus at. Straight ahead. +* target_index: [1] relative index of component to focus at, e.g. next is +1 +* mode_input: [1] order of the phonon mode to be scattered (unphysical). In the range 0 through 11 . If mode_input == 12, then the mode is choosen at random +* verbose_input: [1] toggles verbose mode (verbose>1) or eigenvector mode (verbose==1) +* dispersion: [1] 0: do nothing special; 1: calculate and output the dispersion; 2: calculate and output the dispersion while varying v_f +* +* OUTPUT PARAMETERS: +* V_rho: [AA^-3] Unit cell density +* V_my_s: [m^-1] Attenuation factor due to incoherent scattering +* V_my_a_v: [m^-1] Attenuation factor due to absorbtion +* +* %L +* The test/example instrument Test_Phonon.instr. +* +* %E +******************************************************************************/ + +DEFINE COMPONENT Phonon_BvK_PG +DEFINITION PARAMETERS () +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) +OUTPUT PARAMETERS (q_x, q_y, q_z) +/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ + +SHARE +%{ +//%include "eigen.h" +//%include "eigen.c" +#ifndef PHONON_SIMPLE +#define PHONON_SIMPLE $Revision$ // KL comment 040125: what is this used for? +#pragma acc routine + +// ---------------- THINGS THAT COULD BE GENERAL FOR McStas -------------- + +#define T2E (1/11.605) // Kelvin to meV converter +#define pi 3.14159265358979323846264338327950288419716939937510L +#define sqrt3 1.73205080756887729352744634150587236694280525381038L +#define THz2meV 4.13566 // THz to meV converter + +#define Da2kg 1.66054e-27 //Dalton to kg converter +#define dyn2N 1e-3 // dyn/cm to N/m converter +#define hbar 1.0546E-34 +#define mn (1.0087*Da2kg) + +double nbose(double omega, double T) /* Give it another name ?? */ + { + double nb; + nb= (omega>0) ? 1+1/(exp(omega/(T*T2E))-1) : 1/(exp(-omega/(T*T2E))-1); + return nb; + } +#undef T2E + +/* Routine types from Numerical Recipies book */ +#define UNUSED (-1.11e30) +#define MAXRIDD 60 + +void fatalerror_cpu(char *s) + { + fprintf(stderr,"%s \n",s); + exit(1); + } + + void nrerror(const char *error_text){ + printf("Numerical Recipes run-time error ...\n"); + printf("%s\n",error_text); + printf("... now exiting to system.\n"); + exit(1); +} + +#pragma acc routine +void fatalerror(char *s) + { + #ifndef OPENACC + fatalerror_cpu(s); + #endif + } + + #pragma acc routine + +void bubblesort (int size , double* inputarray , int* index) +//this function modifies the inputarray by bubble sorting, and also modifies the index array as the +//index after the sorting. The index array should be initialized as index[] = {0,1,2,3,4,5,6,7...} +{ + int tempindex = 0; + int i; + int j; + for (i = 0 ; i < size-1 ; i++) + { + for (j = 0 ; j < size-1 ; j++) + { + if (inputarray[index[j]] > inputarray[index[j+1]]) + { + tempindex = index[j]; + index[j] = index[j+1]; + index[j+1] = tempindex; + } + } + } + return; +} + +void mergesort(double *X, int n, int *index) /* Written by James Avery and does not work at present */ +{ + for(int i=0;i=r) continue; /* If right list is empty, sequence is already sorted */ + + int j = 0, j0 = i, j1 = m; + double sorted[segment_length]; + int sorted_index[segment_length]; + + /* Merge the two adjacent sorted lists X[i:m] and X[m:r] */ + double t0 = X[j0], t1 = X[j1]; + while(j +#include +#include +#include +#include "eigen.h" + +extern FILE *matrix_test_file, *parms_test_file; +extern int matrix_test_count; +int mode, verbose, shape; + +#define a_latt 2.461 //PG lattice constant in AA; Nicklow1972 gives 2.45 (WHICH PAGE?) +#define c_latt 6.708 // PG lattice constant in AA; Nicklow1972 gives 6.70 +#define b_length 6.646 // PG scattering legth in fm + +double M = 12.011 * Da2kg; //atomic mass of C + +#define K_l1 3.62e+5 // Force constant longitudinal nn1 - in (a,b) plane; here units of dyn +#define K_t1 1.99e+5 // Force constant transverse nn1 - in (a,b) plane +#define K_l2 1.33e+5 // Force constant longitudinal nn2 - in (a,b) plane +#define K_t2 -0.520e+5 // Force constant transverse nn2 - in (a,b) plane +#define K_l3 -0.037e+5 // Force constant longitudinal nn3 - in (a,b) plane +#define K_t3 0.288e+5 // Force constant transverse nn3 - in (a,b) plane +#define K_l4 0.058e+5 // Force constant longitudinal nn4 - along c +#define K_t4 0.0077e+5 // Force constant transverse nn4 - along c + + //Lattice basis + + // PG is hexagonal close packed with 4 atoms per cell + // Coordinates from Nicklow1972, Fig. 7: Atoms A, B, C, D + double Delta[4*3] = {0 , 0 , 0 , a_latt/(2*sqrt3) , a_latt/2 , 0 , 0 , 0 , c_latt/2 , -a_latt/(2*sqrt3) , a_latt/2 , c_latt/2}; + + // Lattice vectors from paper fig. 7 + double avec[3] = {a_latt*sqrt3/2 , a_latt*0.5 , 0}; // THIS IS NEVER USED ??!! + double bvec[3] = {a_latt*(-sqrt3/2) , a_latt*0.5 , 0}; + double cvec[3] = {0 , 0 , c_latt}; + + // reciprocal lattice vectors + double astar = 4*PI/(sqrt3*a_latt); // length of reciprocal lattice vector a* + double cstar = 2*PI/c_latt; // length of reciprocal lattice vector c* + + // Rotation matrices + // m(12*i+j) = M(i,j) + double Rot120[3][3] = {{-0.5 , sqrt3/2 , 0} , {-sqrt3/2 , -0.5 , 0} , {0 , 0 , 1}}; + double Rot120_flat[9] = {-0.5 , sqrt3/2 , 0 , -sqrt3/2 , -0.5 , 0 , 0 , 0 , 1}; + double Rot60[3][3] = {{0.5 , sqrt3/2 , 0} , {-sqrt3/2 , 0.5 , 0} , {0 , 0 , 1}}; + double Rot120rev[3][3] = {{-0.5 , -sqrt3/2 , 0} , {sqrt3/2 , -0.5 , 0} , {0 , 0 , 1}}; //rotate -120 + double Rot60rev[3][3] = {{0.5 , -sqrt3/2 , 0} , {sqrt3/2 , 0.5 , 0} , {0 , 0 , 1}}; //rotate -60 + double Rot180[3][3] = {{-1, 0, 0}, {0, -1, 0}, {0, 0, 1}}; // rotate 180 or -180 +// double Rot5[3][3] = {{0.5 , sqrt3/2 , 0} , {-sqrt3/2 , 0.5 , 0} , {0 , 0 , 1}}; //rotate 60 +// double Rot5rev[3][3] = {{0.5 , -sqrt3/2 , 0} , {sqrt3/2 , 0.5 , 0} , {0 , 0 , 1}}; //rotate -60 + + //1st neighbour + +// double r_j1[3][3] = {{-a_latt/sqrt3 , 0 , 0} , {a_latt*0.5/sqrt3 , -a_latt*0.5 , 0} , {a_latt*0.5/sqrt3 , a_latt*0.5 , 0}}; // This holds from atoms A and D +// double r_j2[3][3] = {{a_latt/sqrt3 , 0 , 0} , {-a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {-a_latt*0.5/sqrt3 , -a_latt*0.5 , 0}}; // This holds from atoms B and C + + double Phi_nn1[3][3] = {{K_l1*dyn2N , 0 , 0} , {0 , K_t1*dyn2N , 0} , {0 , 0 , K_t1*dyn2N}}; + double Phi1[3][3] = {{K_l1*dyn2N , 0 , 0} , {0 , K_t1*dyn2N , 0} , {0 , 0 , K_t1*dyn2N}}; //Phi1 = Phi_nn1 + double Phi2[3][3]; + double Phi3[3][3]; + + //2nd neighbour +// double r_j11[9][3] = {{a_latt*(-1/sqrt3) , 0 , 0} , {a_latt*0.5/sqrt3 , a_latt*(-0.5) , 0} , {a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {0 , a_latt , 0} , {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , 0}};//r_j1 + r_add +// double r_j21[9][3] = {{-a_latt*(-1/sqrt3) , -a_latt*0 , 0} , {-a_latt*0.5/sqrt3 , -a_latt*(-0.5) , 0} , {-a_latt*0.5/sqrt3 , -a_latt*0.5 , 0} , {0 , a_latt , 0} , {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , 0}};//r_j2 + r_add + + double Phi_nn2[3][3] = {{K_t2*dyn2N , 0 , 0} , {0 , K_l2*dyn2N , 0} , {0 , 0 , K_t2*dyn2N}}; + double Phi4[3][3] = {{K_t2*dyn2N , 0 , 0} , {0 , K_l2*dyn2N , 0} , {0 , 0 , K_t2*dyn2N}}; + double Phi5[3][3]; + double Phi6[3][3]; + double Phi7[3][3]; + double Phi8[3][3]; + double Phi9[3][3]; + + //3rd neighbour +// double r_j12[12][3] = {{a_latt*(-1/sqrt3) , 0 , 0} , {a_latt*0.5/sqrt3 , a_latt*(-0.5) , 0} , {a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {0 , a_latt , 0} , {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , 0} , {2*a_latt/sqrt3 , 0 , 0} , {-a_latt/sqrt3 , a_latt , 0} , {-a_latt/sqrt3 , -a_latt , 0}};//r_j11 + r_add2 + double r_j2[12][3] = {{a_latt/sqrt3 , 0 , 0} , {-a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {-a_latt*0.5/sqrt3 , -a_latt*0.5 , 0} , {0 , a_latt , 0} , {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , 0} , {-2*a_latt/sqrt3 , 0 , 0} , {a_latt/sqrt3 , -a_latt , 0} , {a_latt/sqrt3 , a_latt , 0}};//r_j21 - r_add2 + + double Phi_nn3[3][3] = {{K_l3*dyn2N , 0 , 0} , {0 , K_t3*dyn2N , 0} , {0 , 0 , K_t3*dyn2N}}; + + // c2 = 2/sqrt(6); + // c1 = 1/sqrt(6); + // c0 = 1/sqrt(2); + // c3 = 1/sqrt(3); + + double Phi10[3][3] = {{K_l3*dyn2N , 0 , 0} , {0 , K_t3*dyn2N , 0} , {0 , 0 , K_t3*dyn2N}}; //Phi_nn3 + double Phi11[3][3]; + double Phi12[3][3]; + + //4th neighbour + double r_j13[14][3] = {{-a_latt/sqrt3 , 0 , 0} , {a_latt*0.5/sqrt3 , -a_latt*0.5 , 0} , {a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {0 , a_latt , 0} , {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , 0} , {2*a_latt/sqrt3 , 0 , 0} , {-a_latt/sqrt3 , a_latt , 0} , {-a_latt/sqrt3 , -a_latt , 0} , {0 , 0 , c_latt/2} , {0 , 0 , -c_latt/2}};//r_j12 + r_add3 Sublattice B and D do not have this coupling + + double Phi_nn4[3][3] = {{K_t4*dyn2N , 0 , 0} , {0 , K_t4*dyn2N , 0} , {0 , 0 , K_l4*dyn2N}}; + + double Phi13[3][3] = {{K_t4*dyn2N , 0 , 0} , {0 , K_t4*dyn2N , 0} , {0 , 0 , K_l4*dyn2N}}; //Phi_nn4; + double Phi14[3][3] = {{K_t4*dyn2N , 0 , 0} , {0 , K_t4*dyn2N , 0} , {0 , 0 , K_l4*dyn2N}}; //Phi_nn4; + +// ----------------- END SPECIFIC FOR THE PG COPONENT -------------------------- + + + +/* TODO: Det er meget langsomt at bruge pointere. Erstat med flad memory. */ +// m(12*i+j) = M(i,j) +//3*3 matrices multiplication +void matrixmultiplication( double matrix1[3][3], double matrix2[3][3], double result[3][3]) + { + // multiplies matrix1 by matrix2 and places the result in result +// double matrix1[3][3]={{*(array1), *(array1+1), *(array1+2)}, {*(array1+3), *(array1+4), *(array1+5)}, {*(array1+6), *(array1+7), *(array1+8)}}; +// double matrix2[3][3]={{*(array2), *(array2+1), *(array2+2)}, {*(array2+3), *(array2+4), *(array2+5)}, {*(array2+6), *(array2+7), *(array2+8)}}; + + result[0][0]=matrix2[0][0]*matrix1[0][0]+matrix2[1][0]*matrix1[0][1]+matrix2[2][0]*matrix1[0][2]; + result[0][1]=matrix2[0][1]*matrix1[0][0]+matrix2[1][1]*matrix1[0][1]+matrix2[2][1]*matrix1[0][2]; + result[0][2]=matrix2[0][2]*matrix1[0][0]+matrix2[1][2]*matrix1[0][1]+matrix2[2][2]*matrix1[0][2]; + result[1][0]=matrix2[0][0]*matrix1[1][0]+matrix2[1][0]*matrix1[1][1]+matrix2[2][0]*matrix1[1][2]; + result[1][1]=matrix2[0][1]*matrix1[1][0]+matrix2[1][1]*matrix1[1][1]+matrix2[2][1]*matrix1[1][2]; + result[1][2]=matrix2[0][2]*matrix1[1][0]+matrix2[1][2]*matrix1[1][1]+matrix2[2][2]*matrix1[1][2]; + result[2][0]=matrix2[0][0]*matrix1[2][0]+matrix2[1][0]*matrix1[2][1]+matrix2[2][0]*matrix1[2][2]; + result[2][1]=matrix2[0][1]*matrix1[2][0]+matrix2[1][1]*matrix1[2][1]+matrix2[2][1]*matrix1[2][2]; + result[2][2]=matrix2[0][2]*matrix1[2][0]+matrix2[1][2]*matrix1[2][1]+matrix2[2][2]*matrix1[2][2]; + + return; + } + +// TODO: (James Avery) Benchmark against pointer versions: matrix3mupltiplication or matrixmultiplication +// m(12*i+j) = M(i,j) +void matrix3x3_multiply3(const double A[9], const double B[9], const double C[9], double result[9]) +{ + // R_{i,j} = \sum_{k,l=1}^3 A_{i,k}*B_{k,l}*C_{l,j} + for(int i=0;i<3;i++) + for(int j=0;j<3;j++){ + double ij_sum = 0; + for(int k=0;k<3;k++) + for(int l=0;l<3;l++) + ij_sum += A[i*3+k]*B[k*3+l]*C[l*3+j]; + result[i*3+j] = ij_sum; + } +} + +/* TODO: (James Avery) It is slow to use points. Replace with flat memory. */ + // m(12*i+j) = M(i,j) + // Counter argument: (Kim Lefmann) pointers are now only used in the initialization. In the TRACE code, memory is flat + void matrix3multiplication(double matrix1[3][3], double matrix2[3][3], double matrix3[3][3], double result[3][3]) + { + // Performs the multiplication (matrix1*matrix2)*matrix3 and placed the result in result +//#define TEST_MULTIPLY + double temp[3][3]; + matrixmultiplication(matrix1,matrix2,temp); + matrixmultiplication(temp,matrix3,result); +#ifdef TEST_MULTIPLY + printf("Matrix3multiply called with \n"); + printf("matrix1 = ("); + for (int i=0; i<3; i++) + { + printf("( %g %g %g), ",matrix1[i][0], matrix1[i][1], matrix1[i][2]); + } + printf(")\n Matrix2 = ("); + for (int i=0; i<3; i++) + { + printf("( %g %g %g), ",matrix2[i][0], matrix2[i][1], matrix2[i][2]); + } + printf(")\n Matrix3 = ("); + for (int i=0; i<3; i++) + { + printf("( %g %g %g), ",matrix3[i][0], matrix3[i][1], matrix3[i][2]); + } + printf(")\n Temp = ("); + for (int i=0; i<3; i++) + { + printf("( %g %g %g), ",temp[i][0], temp[i][1], temp[i][2]); + } + printf(")\n Result = ("); + for (int i=0; i<3; i++) + { + printf("( %g %g %g), ",result[i][0], result[i][1], result[i][2]); + } + printf(")\n"); +#endif // TEST_MULTIPLY + + return; + } + + +/* TODO: (James Avery) Replace by standard complex.h operation */ +complex complexexp_qdotr (double* q, double* rj) +{ +// double qrjtest + double qrj = *(q)**(rj)+*(q+1)**(rj+1)+*(q+2)**(rj+2); +// printf(" q= (%g %g %g), r = (%g %g %g) qrj = %g \n",q[0],q[1],q[2],rj[0],rj[1],rj[2],qrj); +// return (cexp(qrj)); // KL comment 030125 WHY does this not run? + return (cos(qrj)+I*sin(qrj)); +} + + +void initialize_omega_q() +{ + // Make the matrix algebra that finalizes the force constant set-up + + matrix3multiplication(Rot120,Phi_nn1,Rot120rev,Phi2); //Phi2=Rot120*Phi_nn1*Rot120^(-1) + matrix3multiplication(Rot120rev,Phi_nn1,Rot120,Phi3); //Phi3=Rot120^{-1}*Phi_nn1*Rot120 + matrix3multiplication(Rot60,Phi_nn2,Rot60rev,Phi5); //Phi5=Rot60*Phi_nn2*Rot60^{-1} + matrix3multiplication(Rot120,Phi_nn2,Rot120rev,Phi6); //Phi6=Rot120*Phi_nn2*Rot120^(-1); + matrix3multiplication(Rot180,Phi_nn2,Rot180,Phi7); //Phi7=Rot180*Phi_nn2*Rot180; + matrix3multiplication(Rot120rev,Phi_nn2,Rot120,Phi8); //Phi8=Rot120^{-1}*Phi_nn2*Rot120; + matrix3multiplication(Rot60rev,Phi_nn2,Rot60,Phi9); //Phi9=Rot60^{-1}*Phi_nn2*Rot60; + matrix3multiplication(Rot120,Phi_nn3,Rot120rev,Phi11);//Phi11=Rot120*Phi_nn3*Rot120^(-1); + matrix3multiplication(Rot120rev,Phi_nn3,Rot120,Phi12);//Phi12=Rot120^{-1}*Phi_nn3*Rot120; + + return; +} + +//---------------------------- START CENTRAL FUNCTION OMEGA-Q ---------------------------- + + double omega_q(complex* parms, complex* eigenvectorgood) { + + // _______________ DETERMINE Q ____________________________ + + double vi, vf, vv_x, vv_y, vv_z, vi_x, vi_y, vi_z; + double q, qx, qy, qz, Jq, hw_phonon, hw_neutron; + double h,k,l,hk_inplane, hk_outofplane; + int disp; + + for(int ii=0;ii<8;ii++) if(isnan(creal(parms[ii])) || isnan(cimag(parms[ii]))){ + fprintf(stderr,"omega_q: Array parms contains a nan at position %d.\n",ii); + fprintf(stderr,"parms = ["); for(int i=0;i<8;i++) fprintf(stderr,"%g + i%g%s",creal(parms[i]), cimag(parms[i]), i+1<8?", ":"]\n"); + abort(); + } + + vf=parms[0]; + vi=parms[1]; + vv_x=parms[2]; + vv_y=parms[3]; + vv_z=parms[4]; + vi_x=parms[5]; + vi_y=parms[6]; + vi_z=parms[7]; + h = parms[9]; + k = parms[10]; + l = parms[11]; + disp = (int)parms[12]; + + qx=V2K*(vi_x-vf*vv_x); + qy=V2K*(vi_y-vf*vv_y); + qz=V2K*(vi_z-vf*vv_z); + + if (disp==1) /* calculate dispersion directly */ + { + printf("Dispersion calculation: "); + qx = h*astar; + qy = k*astar; //correct only for qy=0 ! + qz = l*cstar; + } + + double qvec[3]={qx,qy,qz}; + + if (verbose>=4) + printf("Enter omega_q; q = (%g, %g, %g) \n",qx,qy,qz); + + // ________________ END DETERMINE Q ________________- + + + complex Phi_diag[3*3]; + int i; + int j; + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + Phi_diag[i*3+j] = Phi1[i][j]+Phi2[i][j]+Phi3[i][j]+Phi4[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[3][0]))+Phi5[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[4][0]))+Phi6[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[5][0]))+Phi7[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[6][0]))+Phi8[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[7][0]))+Phi9[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[8][0]))+Phi10[i][j]+Phi11[i][j]+Phi12[i][j]; + } + } + + complex Phi_diag1[3*3]; + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + Phi_diag1[i*3+j] = Phi13[i][j]+Phi14[i][j]; + } + } + + complex Phi_offdiag1[3*3]; + complex Phi_offdiag2[3*3]; + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + Phi_offdiag1[i*3+j] =-Phi1[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[0][0]))-Phi2[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[1][0]))-Phi3[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[2][0]))-Phi10[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[9][0]))-Phi11[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[10][0]))-Phi12[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[11][0])); + Phi_offdiag2[i*3+j] = conj(Phi_offdiag1[i*3+j]); + } + } + + complex Phi_6Dim[6*6]; + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + Phi_6Dim[i*6+j] = Phi_diag[i*3+j]+Phi_diag1[i*3+j]; + Phi_6Dim[i*6+j+3] = Phi_offdiag1[i*3+j]; + Phi_6Dim[(i+3)*6+j] = Phi_offdiag2[i*3+j]; + Phi_6Dim[(i+3)*6+j+3] = Phi_diag[i*3+j]; // Because Phi_diag1 does not appear for the B,C, sublattices + } + } + + complex Phi_AC[3*3]; + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + Phi_AC[i*3+j] = -Phi13[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[12][0]))-Phi14[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[13][0])); + } + } + +// complex Zero_3Dim[3][3] = {{0 , 0 , 0},{0 , 0 , 0},{0 , 0 , 0}}; + + complex Phi_6Dim_offdiag[6*6]; + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + Phi_6Dim_offdiag[i*6+j] = Phi_AC[i*3+j]; + Phi_6Dim_offdiag[i*6+j+3] = (complex)0.0; + Phi_6Dim_offdiag[(i+3)*6+j] = (complex)0.0; + Phi_6Dim_offdiag[(i+3)*6+j+3] = (complex)0.0; + } + } + + complex Phi_12Dim[12*12]; + for (i = 0; i < 6; i++) { + for (j = 0; j < 6; j++) { + Phi_12Dim[i*DIM+j] = Phi_6Dim[i*6+j]; + Phi_12Dim[i*DIM+j+6] = Phi_6Dim_offdiag[i*6+j]; + Phi_12Dim[(i+6)*DIM+j] = conj(Phi_6Dim_offdiag[i*6+j]); + Phi_12Dim[(i+6)*DIM+j+6] = Phi_6Dim[i*6+j]; + } + } + +#ifdef TEST_MATRICES + printf("Phi_12Dim=\n"); + for (i = 0; i < 12; i++) { + for (j = 0; j< 12; j++) { + if (j == 0) { + printf("{"); + } + printf("%g+(%g)i ", creal(Phi_12Dim[i*DIM+j]) , cimag((Phi_12Dim[i*DIM+j]))); + if (j == 11) { + printf("}\n"); + } + } + } + + printf("Phi_6Dim=\n"); + for (i = 0; i < 6; i++) { + for (j = 0; j< 6; j++) { + if (j == 0) { + printf("{"); + } + printf("%g+(%g)i ", creal(Phi_6Dim[i*6+j]) , cimag((Phi_6Dim[i*6+j]))); + if (j == 5) { + printf("}\n"); + } + } + } + + printf("Phi_diag=\n"); + for (i = 0; i < 3; i++) { + for (j = 0; j< 3; j++) { + if (j == 0) { + printf("{"); + } + printf("%g+(%g)i ", creal(Phi_diag[i*3+j]) , cimag((Phi_diag[i*3+j]))); + if (j == 2) { + printf("}\n"); + } + } + } + printf("Phi2={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}}\n",Phi2[0][0],Phi2[0][1],Phi2[0][2],Phi2[1][0],Phi2[1][1],Phi2[1][2],Phi2[2][0],Phi2[2][1],Phi2[2][2]); + printf("Phi3={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi3[0][0],Phi3[0][1],Phi3[0][2],Phi3[1][0],Phi3[1][1],Phi3[1][2],Phi3[2][0],Phi3[2][1],Phi3[2][2]); + printf("Phi5={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi5[0][0],Phi5[0][1],Phi5[0][2],Phi5[1][0],Phi5[1][1],Phi5[1][2],Phi5[2][0],Phi5[2][1],Phi5[2][2]); + printf("Phi6={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi6[0][0],Phi6[0][1],Phi6[0][2],Phi6[1][0],Phi6[1][1],Phi6[1][2],Phi6[2][0],Phi6[2][1],Phi6[2][2]); + printf("Phi7={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi7[0][0],Phi7[0][1],Phi7[0][2],Phi7[1][0],Phi7[1][1],Phi7[1][2],Phi7[2][0],Phi7[2][1],Phi7[2][2]); + printf("Phi8={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi8[0][0],Phi8[0][1],Phi8[0][2],Phi8[1][0],Phi8[1][1],Phi8[1][2],Phi8[2][0],Phi8[2][1],Phi8[2][2]); + printf("Phi9={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi9[0][0],Phi9[0][1],Phi9[0][2],Phi9[1][0],Phi9[1][1],Phi9[1][2],Phi9[2][0],Phi9[2][1],Phi9[2][2]); + printf("Phi11={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi11[0][0],Phi11[0][1],Phi11[0][2],Phi11[1][0],Phi11[1][1],Phi11[1][2],Phi11[2][0],Phi11[2][1],Phi11[2][2]); + printf("Phi12={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi12[0][0],Phi12[0][1],Phi12[0][2],Phi12[1][0],Phi12[1][1],Phi12[1][2],Phi12[2][0],Phi12[2][1],Phi12[2][2]); + printf("Phi14={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi14[0][0],Phi14[0][1],Phi14[0][2],Phi14[1][0],Phi14[1][1],Phi14[1][2],Phi14[2][0],Phi14[2][1],Phi14[2][2]); + + printf(" 1-exp(q.r_j13[3]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[3][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[3][0]))); + printf(" 1-exp(q.r_j13[4]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[4][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[4][0]))); + printf(" 1-exp(q.r_j13[5]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[5][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[5][0]))); + printf(" 1-exp(q.r_j13[6]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[6][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[6][0]))); + printf(" 1-exp(q.r_j13[7]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[7][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[7][0]))); + printf(" 1-exp(q.r_j13[8]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[8][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[8][0]))); + +#endif // TEST_MATRICES + + double eigenvalue[DIM]; + complex Q[DIM*DIM]; + complex matrix[DIM*DIM]; + int index[DIM]; + + + for (i=0; i= 6) + { + printf("Before bubblesort\n"); + for (i=0; i= 5) + { + printf("After bubblesort\n"); + for (i=0; i=7) + { + printf("Q = ["); + for (i = 0; i < DIM ; i++) + { + printf("\n ("); + for (j=0; j=5) + { + printf("Diagonalization done, eigenvalues (energies) are: \n ("); + for (i=0; i= 4){ + printf("Return from omega_q \n"); + printf("hw_neutron = %g (vi = %g, vf = %g)\n",hw_neutron, vi, vf); + printf("hw_phonon = %g = eigenvaluegood[%d] = %g = eigenvalues[%d] = %g\n", + hw_phonon,mode,eigenvaluegood[mode],index[mode], eigenvalue[index[mode]]); +// printf("old parms = ["); for(int i=0;i<8;i++) fprintf(stderr,"%g+i%g%s",creal(parms[i]), cimag(parms[i]), i+1<8?", ":"]\n"); + } + + parms[8] = hw_phonon; + + if (disp==1) { + printf("mode = %d, (h,k,l) = (%g,%g,%g), hw_q = %g",mode,h,k,l,hw_phonon); + printf(" polarization: ( "); + for (j = 0; j < DIM; j++) + { + printf("%g + i %g ,",creal(eigenvectorgood[j]),cimag(eigenvectorgood[j])); + } + printf(")\n"); + } + if (disp==2) { + hk_inplane = qx/astar; + hk_outofplane = qy/astar; // Derivation assume astar along x and 60 degrees between astar and bstar: + // h_ip = h + k/2 h_op = sqrt(3)*k/2 => k = 2/sqrt(3)*h_op h = h_ip - h_op/sqrt(3) + k = 2/sqrt(3)*hk_outofplane; + h = hk_inplane - hk_outofplane/sqrt(3); + l = qz/cstar; + printf("mode = %d, vf = %g, (h,k,l) = (%g,%g,%g), hw_q = %g",mode,vf,h,k,l,hw_phonon); + printf(" polarization: (%g + i %g , %g + i %g , %g + i %g) \n",creal(eigenvectorgood[0]),cimag(eigenvectorgood[0]),creal(eigenvectorgood[1]),cimag(eigenvectorgood[1]), creal(eigenvectorgood[2]),cimag(eigenvectorgood[2])); + } + + return (hw_phonon - hw_neutron); +} + +//--------------------- END CENTRAL FUNCTION OMEGA-Q ------------------------ + +// ------------------START OF THE GENERAL ROOT-FINDING CODE-------------------- + +double last_fh, last_x2; +double zridd(double (*func)(complex*,complex*), double x1, double x2, complex *parms, complex* vector, double xacc) + { + int j; + double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; + + // fprintf(stderr,"i zridd(%g,%g, parms,%g)\n",x1,x2,xacc); + parms[0]=x1; + // fprintf(stderr,"zridd fl = omega_q\n"); + if (xacc>0) + { + fl=(*func)(parms,vector); + } + else + { + fl = last_fh; + xacc = -xacc; + if (fabs(x1-last_x2)>xacc) // KL comment 030125 WHY does the library routine zridd use complex abs() function ?? changed to fabs() + printf("*** error in zridd *** x1: %g last_x2: %g \n",x1,last_x2); + } + // fprintf(stderr,"fl = %g\n",fl); + parms[0]=x2; + // fprintf(stderr,"zridd fh = omega_q; parms[0] = x2 = %g+i%g = %g\n",creal(parms[0]), cimag(parms[0]),x2); + last_fh=fh=(*func)(parms,vector); + last_x2 = x2; + // fprintf(stderr,"fh = %g\n",fh); + if (fl*fh >= 0) + { + if (fl==0) return x1; + if (fh==0) return x2; + return UNUSED; + } + else + { + xl=x1; + xh=x2; +// printf("zridd sign change: v_low : %g v_high: %g \n",xl,xh); + ans=UNUSED; + for (j=1; j= fh ? 1.0 : -1.0)*fm/s); + if (fabs(xnew-ans) <= xacc) + return ans; + ans=xnew; + parms[0]=ans; + // fprintf(stderr,"s = %g; fl, fm, fh = %g,%g,%g; fnew = %g; fm*fm-fl*fh = %g\n", + // s, fl, fm,fh, fnew, fm*fm-fl*fh); + //fprintf(stderr,"zridd fnew = omega_q\n"); + fnew=(*func)(parms,vector); + if (fnew == 0.0) return ans; + if (fabs(fm)*SIGN(fnew) != fm) + { + xl=xm; + fl=fm; + xh=ans; + fh=fnew; + } + else + if (fabs(fl)*SIGN(fnew) != fl) + { + xh=ans; + fh=fnew; + } + else + if(fabs(fh)*SIGN(fnew) != fh) + { + xl=ans; + fl=fnew; + } + else + fatalerror("never get here in zridd"); + if (fabs(xh-xl) <= xacc) + return ans; + } + fatalerror("zridd exceeded maximum iterations"); + } + return 0.0; /* Never get here */ + } + +#pragma acc routine +double zridd_gpu(double x1, double x2, complex *parms, complex* vector, double xacc) + { + int j; + double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; + + parms[0]=x1; + // fprintf(stderr,"fl=omega_q(%g+i%g)\n",creal(parms[0]),cimag(parms[0])); + fl=omega_q(parms,vector); + parms[0]=x2; + // fprintf(stderr,"fh=omega_q(%g+i%g)\n",creal(parms[0]),cimag(parms[0])); + fh=omega_q(parms,vector); + if (fl*fh >= 0) + { + if (fl==0) return x1; + if (fh==0) return x2; + return UNUSED; + } + else + { + xl=x1; + xh=x2; + ans=UNUSED; + for (j=1; j= fh ? 1.0 : -1.0)*fm/s); + if (fabs(xnew-ans) <= xacc) + return ans; + ans=xnew; + parms[0]=ans; + // fprintf(stderr,"fnew=omega_q(%g+i%g)\n",creal(parms[0]),cimag(parms[0])); + fnew=omega_q(parms,vector); + if (fnew == 0.0) return ans; + if (fabs(fm)*SIGN(fnew) != fm) + { + xl=xm; + fl=fm; + xh=ans; + fh=fnew; + } + else + if (fabs(fl)*SIGN(fnew) != fl) + { + xh=ans; + fh=fnew; + } + else + if(fabs(fh)*SIGN(fnew) != fh) + { + xl=ans; + fl=fnew; + } + else + fatalerror("never get here in zridd"); + if (fabs(xh-xl) <= xacc) + return ans; + } + fatalerror("zridd exceeded maximum iterations"); + } + return 0.0; /* Never get here */ + } + +#define ROOTACC 1e-8 + + int findroots(double brack_low, double brack_mid, double brack_high, double *list, int* index, double (*f)(complex*,complex*), complex *parms, complex* vector, int low_steps, int high_steps) + // *f is here omega_q + // *parms is here p_call + { + double root,range; + int i; + + if (verbose==2) + printf("Enter findroots() \n"); + + // First, find roots in energy loss side, if any + range = brack_mid-brack_low; + for (i=0; i<(low_steps-1); i++) + { + if (i==0) + root = zridd(f, brack_low+range*i/low_steps, + brack_low+range*(i+1)/low_steps, + (complex *)parms, (complex *)vector, ROOTACC); + else + root = zridd(f, brack_low+range*i/low_steps, + brack_low+range*(i+1)/low_steps, + (complex *)parms, (complex *)vector, -ROOTACC); // Re-use the previous fh value + if (root != UNUSED) + { + list[(*index)++]=root; + if (verbose >=4) + printf("--- findroots returned a root on the energy loss side: vf = %g \n",root); + } + } + + range = (brack_mid-brack_low)/(double)low_steps; + for (i=0; i=4) + printf("--- findroots returned a root on the energy loss side: vf = %g \n",root); + } + } + + // Second, find roots in energy gain side, there is always some + range = (brack_high-brack_mid)/(double)high_steps; + for (i=0; i= 4) + printf("*** findroots returned a root on the energy gain side: vf = %g \n",root); + } + } + + range = brack_high-brack_mid; + for (i=1; i= 4) + printf("*** findroots returned a root on the energy gain side: vf = %g \n",root); + } + } + + } + +#pragma acc routine + int findroots_gpu(double brack_low, double brack_mid, double brack_high, double *list, int* index, double *parms, complex* vector, int low_steps, int high_steps) + { + double root,range=brack_mid-brack_low; + int i; + + for (i=0; i 0 && yheight) shape=0; /* cylinder */ + + + V_rho = 1/(a_latt*a_latt*c_latt*sqrt3/2); // (unit: Å^-3) + V_my_s = (V_rho * 100 * sigma_inc); + V_my_a_v = (V_rho * 100 * sigma_abs * 2200); + + /* now compute target coords if a component index is supplied */ + if (!target_index && !target_x && !target_y && !target_z) target_index=1; + if (target_index){ + Coords ToTarget; + ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+target_index), POS_A_CURRENT_COMP); + ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget); + coords_get(ToTarget, &target_x, &target_y, &target_z); + } + if (!(target_x || target_y || target_z)) { + printf("Phonon_BkV_PG: %s: The target is not defined. Using direct beam (Z-axis).\n", + NAME_CURRENT_COMP); + target_z=1; + } + initialize_omega_q(); + verbose = verbose_input; + + // m(12*i+j) = M(i,j) + // OK for: Rot120 + for (i=0; i<3; i++) + for (j=0; j<3; j++) + { + printf("i,j: %i, %i ROT120: %g ROT120FLAT; %g \n",i,j,Rot120[i][j],Rot120_flat[3*i+j]); + } + %} + +TRACE +%{ + double t0, t1, t2, t3; /* Entry/exit times for shape */ + double v_i, v_f; /* Neutron velocities: initial, final */ + double vx_i, vy_i, vz_i; /* Neutron initial velocity vector */ + double dt0, dt; /* Flight times through sample */ + double l_full; /* Flight path length for non-scattered neutron */ + double l_i, l_o; /* Flight path lenght in/out for scattered neutron */ + double my_a_i; /* Initial attenuation factor */ + double my_a_f; /* Final attenuation factor */ + double solid_angle; /* Solid angle of target as seen from scattering point */ + double aim_x=0, aim_y=0, aim_z=1; /* Position of target relative to scattering point */ + double omega; /* Energy transfer */ + double qsquare; /* Square of the scattering vector */ + double bose_factor; /* Calculated value of the Bose factor */ + int nf, index; /* Number of allowed final velocities */ + double vf_list[20]; /* List of allowed final velocities */ + double J_factor, J0; /* Jacobian from delta fnc.s in cross section; elastic J-factor */ + double f1, f2; /* Probed values of omega_q minus omega */ + double p_att,p_MC,p_ph1,p_ph2,p_ph3,p_J, p_tot; /* Temporary multipliers */ + complex Gn; /* Phonon structure factor */ + complex q_dot_e; /* q dot e */ + double q_dot_r; /* q dot r */ + double qh,qk,ql,qhkx,qhky,qh_Bragg,qk_Bragg,ql_Bragg,qh_res,qk_res,qh_add,qk_add; /* components of q */ + double dist_00, dist_01, dist_10, dist_11, dist_min; /* used to calculate nearest Bragg point */ + int i,j, intersect; + + if (shape == 0) + intersect = cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius, yheight); + else if (shape == 1) + intersect = box_intersect(&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth); + + if(intersect) + { + // fprintf(stderr,"intersect\n"); + if(t0 < 0) + ABSORB; /* Neutron came from the sample or begins inside */ + + if (verbose>=2) + printf("neutron entered Phonon_BvK_PG \n"); + /* Neutron enters at t=t0. */ + dt0 = t3-t0; /* Time in sample */ + v_i = sqrt(vx*vx + vy*vy + vz*vz); + l_full = v_i * dt0; /* Length of path through sample if not scattered */ + dt = rand01()*dt0; /* Time of scattering (relative to t0) */ + l_i = v_i*dt; /* Penetration in sample at scattering */ + vx_i=vx; + vy_i=vy; + vz_i=vz; + PROP_DT(dt+t0); /* Point of scattering */ + + aim_x = target_x-x; /* Vector pointing at target (e.g. analyzer) */ + aim_y = target_y-y; + aim_z = target_z-z; + + if(focus_aw && focus_ah) { + randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle, + aim_x, aim_y, aim_z, focus_aw, focus_ah, ROT_A_CURRENT_COMP); + } else if(focus_xw && focus_yh) { + randvec_target_rect(&vx, &vy, &vz, &solid_angle, + aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP); + } else { + randvec_target_sphere(&vx,&vy,&vz,&solid_angle,aim_x,aim_y,aim_z, focus_r); + } + NORM(vx, vy, vz); + nf=0; + if (mode_input == 12) + mode = round(12*rand01()-0.5); // Select the phonon mode + else + mode = mode_input; // Use the given mode + if ((mode < 0) || (mode > 11)) // Undefined mode specification + { + printf("mode = %d ",mode); + nrerror("illegal value of mode"); + } + + p_call[0]=-1; // in the end this will be v_f + p_call[1]=v_i; + p_call[2]=vx; + p_call[3]=vy; + p_call[4]=vz; + p_call[5]=vx_i; + p_call[6]=vy_i; + p_call[7]=vz_i; + p_call[9]=hh; + p_call[10]=kk; + p_call[11]=ll; + p_call[12] = (double)dispersion; + + if (dispersion==1) + { + f1=omega_q(p_call,eigenvectormode); + p_call[12]=0.0; + } + +/* if (dispersion==2) + { + for(i=0; i=2) + printf("Call findroots \n"); + findroots(0, v_i, v_i+V_HIGH, vf_list, &nf, omega_q, p_call, eigenvectormode, e_steps_low, e_steps_high); + if (verbose >=2) + printf( "Findroots returned %d roots \n",nf); + #else + findroots_gpu(0, v_i, v_i+V_HIGH, vf_list, &nf, p_call, eigenvectormode, e_steps_low, e_steps_high); + #endif + + if (verbose>=3) + { + printf("Findroots done, mode %i , last phonon energy %g \n",mode,creal(p_call[8])); + } + + // ________________ ALL FROM HERE IS CALCULATION OF INTENSITY - SHOULD BE REVISITED ___________________ + + index=(int)floor(rand01()*nf); // Select random solution + v_f=vf_list[index]; + p_call[0]=v_f; // transfer choice of v_f to omega_q + + f1=omega_q(p_call,eigenvectormode); + if (verbose >= 2) + { + printf("Chosen solution is %i; v_f = %g, hw = %g, energy match is %g; eigenvector is: (", index, v_f, creal(p_call[8]), f1); + for (j=0; j= 3) + { + printf("q transforms in r.l.u to (h, k, l) = (%g, %g, %g). Nearest Bragg point: (%g, %g, %g)\n",qh,qk,ql,qh_Bragg,qk_Bragg,ql_Bragg); + printf("--------------------------------------------\n"); + } + + + if (dispersion==1) + printf("Dispersion found: (h,k,l) = (%g, %g, %g), hw = %g \n",q[0]/astar,q[1]/astar,q[2]/cstar,omega); + + qsquare=1; + Gn = 0; + + double q_Bragg[3]; // Bragg point in Å-1; Cartesian coordinates + q_Bragg[0]=qh_Bragg*astar+qk_Bragg*astar/2.0; + q_Bragg[1]=qk_Bragg*astar*2.0/sqrt3; + q_Bragg[2]=ql_Bragg*cstar; + + for (i=0; i<4; i++) // Calculate phonon structure factor; loop over atoms in unit cell + { + q_dot_e = (complex)0; + q_dot_r = 0; + for (j=0; j<3; j++) // loop over x,y,z + { + q_dot_e += (complex)q[j]*eigenvectormode[3*i+j]; + q_dot_r += q_Bragg[j]*Delta[3*i+j]; + } + if (verbose >= 7) + { + printf("i = %i , q dot r = %g , q dot e = (%g + i %g) Fn contrib = (%g + i %g)",i,q_dot_r,creal(q_dot_e),cimag(q_dot_e),creal(b_length*q_dot_e*cexp(I*q_dot_r)),cimag(b_length*q_dot_e*cexp(I*q_dot_r))); + } + Gn += b_length*q_dot_e*cexp(I*q_dot_r); // Unit fm Å-1 ; for the general lattice, this should be divided by sqrt(mass[j]) + } + if (verbose >= 7) + printf("\n"); + + if (shape == 0) + intersect = cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius, yheight); + else if (shape == 1) + intersect = box_intersect(&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth); + + if(!intersect) + { + printf("FATAL ERROR: Did not hit sample surface from inside.\n"); + exit(1); + } + + if (verbose>=2) + printf("neutron left Phonon_BvK_PG, p_init = %g,",p); + + dt = t3; + l_o = v_f*dt; + my_a_i = V_my_a_v/v_i; + my_a_f = V_my_a_v/v_f; + bose_factor=nbose(omega,T); + J0 = hbar*hbar*v_f*V2K*1E10/mn; // The J-factor for a flat mode + + + p_att = exp(-(V_my_s*(l_i+l_o)+my_a_i*l_i+my_a_f*l_o)); /* Absorption factor (units: 1) */ + p_MC = nf*solid_angle*l_full; /* Focusing factors; assume random choice of n_f possibilities (units: m) */ + if (mode_input == 12) + p_MC *= 12; //The factor 12 is due to MC choice between the 12 modes + p_ph1 = (v_f/v_i)*DW*bose_factor*V_rho*1E30; /* Cross section factor 1 (units: m^-3) */ + p_ph2 = hbar*hbar/(2*(fabs(omega)*1.6022E-22)); /* Jacobian of delta functions in cross section; and constants. Units: J s^2 */ + p_ph3 = (Gn*conj(Gn))*1E-10/M; /* Structure factor. (units: kg^-1) */ + p_J = J0/J_factor; + p_tot = p_att*p_MC*p_ph1*p_ph2*p_ph3*p_J; + p *= p_tot; + + if (verbose>=2) { + printf("nf = %d, solid_angle = %g, l_full = %g \n",nf,solid_angle,l_full); + printf("(p in SI units): p_att = %g, p_MC= %g, p_ph1 = %g, p_ph2 = %g, p_ph3 = %g, p_J = %g, bose_factor = %g, V_rho = %g, J_factor/J0 = %g, G_factor = %g +i %g, conj(G) = %g + I %g, |G|^2 = %g + I %g, ",p_att, p_MC, p_ph1, p_ph2, p_ph3, p_J, bose_factor, V_rho, J_factor/J0, creal(Gn), cimag(Gn), creal(conj(Gn)), cimag(conj(Gn)), creal(Gn*conj(Gn)), cimag(Gn*conj(Gn)) ); + printf(" p_tot = %g, p_final = %g\n",p_tot, p) ; + printf("q is (h,k,l) = (%g, %g, %g) \n",q_x/astar,q_y/astar,q_z/cstar); + } + } /* else transmit: Neutron did not hit the sample */ +%} + +MCDISPLAY +%{ + if (radius > 0 && yheight) { /* cylinder */ + cylinder(0,0,0,radius,yheight,0, 0, 1, 0); + } + 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 diff --git a/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/Test_Phonon_BvK_PG.instr b/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/Test_Phonon_BvK_PG.instr new file mode 100644 index 000000000..187dcef05 --- /dev/null +++ b/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: +* +* %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 +* +* +* %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 diff --git a/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/eigen.c b/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/eigen.c new file mode 100644 index 000000000..03504de7f --- /dev/null +++ b/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/eigen.c @@ -0,0 +1,451 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +/* TODO: Tridiagonal input/outpout in (3,m) storage. */ +/* TODO: Symmetric tridiagonal input/outpout in (2,m) storage. */ +/* TODO: Add pivot to handle ill-conditioned matrices */ +#include "eigen.h" + +#define G PRINTF_G +#define ROWS 0 +#define COLS 1 + +/* TOOD: Generic matrix functions to separate file */ +INLINE real_t vector_norm2(const scalar *x, int n){ + real_t sum = 0; + for(int i=0;i numerical_zero) // Only process if subdiagonal element is not already zero. + { + a[0] = D[k]; a[1] = L[k]; // a = T[k:k+2,k] is the vector of nonzeros in kth subdiagonal column. + + real_t anorm = sqrt(a[0]*a[0] + a[1]*a[1]); + + // // Udrullet + // // reflection_vector(a,anorm,v); + v[0] = D[k]; v[1] = L[k]; + real_t alpha = -copysign(anorm,a[0]); // Koster ingenting + v[0] -= alpha; + + real_t vnorm = sqrt(v[0]*v[0]+v[1]*v[1]); + real_t norm_inv = 1/vnorm; /* Normalize */ + v[0] *= norm_inv; v[1] *= norm_inv; + + Vout[2*k] = v[0]; Vout[2*k+1] = v[1]; + + // // Udrullet + // // apply_reflection(T({k,k+2},{k,k+3}),v); + // // if(k+1=0;k--){ + // We start by targeting the (n,n)-eigenvalue, and gradually + // deflate, working on smaller and smaller submatrices. + real_t d = D[k]; // d = T(k,k) + real_t shift = d; + + // The Gershgorin disk radius is defined by just the row-sums of + // absolute off-diagonal elements, since T is symmetric. As T is + // tridiagonal, only T(k,k-1),T(k,k), and T(k,k+1) are nonzero. + // Thus, the k'th Gershgorin radius is just |T(k,k-1)| + + // |T(k,k+1)| = |T(k,k-1)| + |T(k+1,k)| = |L[k-1]|+|L[k]|. + int i=0; + + real_t GR = (k>0?fabs(L[k-1]):0)+fabs(L[k]); // Initial Gershgorin radius + int not_done = 2; + while(not_done > 0){ + i++; + T_QTQ(k+1, D,L, D,L, V, shift); + if(Q!=0) apply_real_reflections(V,k,Qhat,n); + + GR = (k>0?fabs(L[k-1]):0)+(k+10){ + real_pair l = eigvalsh2x2(D[k-1],L[k-1], /* T[(k-1):k, (k-1):k] */ + L[k-1],D[k] ); + real_t l0 = l.value[0], l1 = l.value[1]; + + shift = fabs(l0-d) < fabs(l1-d)? l0 : l1; // Pick closest eigenvalue + } else + shift = D[k]; + + if(GR <= gershgorin_tolerance) not_done--; + if(i>max_iterations && GR > gershgorin_tolerance){ + fprintf(stderr,"Cannot converge eigenvalue %d to tolerance " G + " in %d iterations using machine precision " G " (d=" G ", shift=" G ", GR=" G ")\n" + "D[k] = " G ", L[k-1] = " G ", L[k] = " G "\n", + k,gershgorin_tolerance,i, + machine_precision,d,shift,GR, + D[k], (k>0)?L[k-1]:0, (k+1 +#include +#include +#include + +#define min(a,b) ((a)<(b)?(a):(b)) + +void mergesort(double *X, int n, int *index) +{ + for(int i=0;i=r) continue; /* If right list is empty, sequence is already sorted */ + +#if TEST_SORT + printf("Step %d:%d: Merging pre-sorted segments: [", step,i); + for(int k=i;k +#include + +#define FLOAT_TYPE float64 + +#if FLOAT_TYPE==float64 + #define real_t double + #define machine_precision DBL_EPSILON + #define PRINTF_G "%g" + +#elif FLOAT_TYPE==float80 + #define real_t long double + #define PRINTF_G "%Lg" + #define creal(x) creall(x) + #define cimag(x) cimagl(x) + + #define conj conjl + #define sqrt sqrtl + #define fabs fabsl + #define machine_precision LDBL_EPSILON + +#elif FLOAT_TYPE==float32 + #define real_t float + #define PRINTF_G "%g" + #define creal(x) crealf(x) + #define cimag(x) cimagf(x) + + #define conj conjf + #define sqrt sqrtf + #define fabs fabsf + #define machine_precision FLT_EPSILON +#endif + +#ifdef __cplusplus + #include + typedef std::complex complex_t; +#else +typedef real_t complex complex_t; +#endif + + +#define INLINE inline __attribute__((always_inline)) + + +typedef struct { + real_t value[2]; +} real_pair; + +typedef complex_t scalar; From fd9ff3babf209ff3742b1ff2fc16551b043c3a4c Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Tue, 24 Feb 2026 13:30:26 +0100 Subject: [PATCH 02/15] Suppress (unused) Avery mergesort in Phonon_BvK_PG. (Mismatching prototype vs. system-provided sort on macOS) --- mcstas-comps/contrib/Phonon_BvK_PG.comp | 40 ------------------------- 1 file changed, 40 deletions(-) diff --git a/mcstas-comps/contrib/Phonon_BvK_PG.comp b/mcstas-comps/contrib/Phonon_BvK_PG.comp index c2e93cdb4..bb2812190 100644 --- a/mcstas-comps/contrib/Phonon_BvK_PG.comp +++ b/mcstas-comps/contrib/Phonon_BvK_PG.comp @@ -150,46 +150,6 @@ void bubblesort (int size , double* inputarray , int* index) return; } -void mergesort(double *X, int n, int *index) /* Written by James Avery and does not work at present */ -{ - for(int i=0;i=r) continue; /* If right list is empty, sequence is already sorted */ - - int j = 0, j0 = i, j1 = m; - double sorted[segment_length]; - int sorted_index[segment_length]; - - /* Merge the two adjacent sorted lists X[i:m] and X[m:r] */ - double t0 = X[j0], t1 = X[j1]; - while(j Date: Tue, 24 Feb 2026 13:43:56 +0100 Subject: [PATCH 03/15] Removed unused mergesort from J Avery --- .../Tests_samples/Test_Phonon_BvK_PG/eigen.c | 451 ------------------ .../Tests_samples/Test_Phonon_BvK_PG/eigen.h | 40 -- .../Test_Phonon_BvK_PG/mergesort.c | 73 --- .../Tests_samples/Test_Phonon_BvK_PG/types.h | 50 -- 4 files changed, 614 deletions(-) delete mode 100644 mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/eigen.c delete mode 100644 mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/eigen.h delete mode 100644 mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/mergesort.c delete mode 100644 mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/types.h diff --git a/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/eigen.c b/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/eigen.c deleted file mode 100644 index 03504de7f..000000000 --- a/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/eigen.c +++ /dev/null @@ -1,451 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include - -/* TODO: Tridiagonal input/outpout in (3,m) storage. */ -/* TODO: Symmetric tridiagonal input/outpout in (2,m) storage. */ -/* TODO: Add pivot to handle ill-conditioned matrices */ -#include "eigen.h" - -#define G PRINTF_G -#define ROWS 0 -#define COLS 1 - -/* TOOD: Generic matrix functions to separate file */ -INLINE real_t vector_norm2(const scalar *x, int n){ - real_t sum = 0; - for(int i=0;i numerical_zero) // Only process if subdiagonal element is not already zero. - { - a[0] = D[k]; a[1] = L[k]; // a = T[k:k+2,k] is the vector of nonzeros in kth subdiagonal column. - - real_t anorm = sqrt(a[0]*a[0] + a[1]*a[1]); - - // // Udrullet - // // reflection_vector(a,anorm,v); - v[0] = D[k]; v[1] = L[k]; - real_t alpha = -copysign(anorm,a[0]); // Koster ingenting - v[0] -= alpha; - - real_t vnorm = sqrt(v[0]*v[0]+v[1]*v[1]); - real_t norm_inv = 1/vnorm; /* Normalize */ - v[0] *= norm_inv; v[1] *= norm_inv; - - Vout[2*k] = v[0]; Vout[2*k+1] = v[1]; - - // // Udrullet - // // apply_reflection(T({k,k+2},{k,k+3}),v); - // // if(k+1=0;k--){ - // We start by targeting the (n,n)-eigenvalue, and gradually - // deflate, working on smaller and smaller submatrices. - real_t d = D[k]; // d = T(k,k) - real_t shift = d; - - // The Gershgorin disk radius is defined by just the row-sums of - // absolute off-diagonal elements, since T is symmetric. As T is - // tridiagonal, only T(k,k-1),T(k,k), and T(k,k+1) are nonzero. - // Thus, the k'th Gershgorin radius is just |T(k,k-1)| + - // |T(k,k+1)| = |T(k,k-1)| + |T(k+1,k)| = |L[k-1]|+|L[k]|. - int i=0; - - real_t GR = (k>0?fabs(L[k-1]):0)+fabs(L[k]); // Initial Gershgorin radius - int not_done = 2; - while(not_done > 0){ - i++; - T_QTQ(k+1, D,L, D,L, V, shift); - if(Q!=0) apply_real_reflections(V,k,Qhat,n); - - GR = (k>0?fabs(L[k-1]):0)+(k+10){ - real_pair l = eigvalsh2x2(D[k-1],L[k-1], /* T[(k-1):k, (k-1):k] */ - L[k-1],D[k] ); - real_t l0 = l.value[0], l1 = l.value[1]; - - shift = fabs(l0-d) < fabs(l1-d)? l0 : l1; // Pick closest eigenvalue - } else - shift = D[k]; - - if(GR <= gershgorin_tolerance) not_done--; - if(i>max_iterations && GR > gershgorin_tolerance){ - fprintf(stderr,"Cannot converge eigenvalue %d to tolerance " G - " in %d iterations using machine precision " G " (d=" G ", shift=" G ", GR=" G ")\n" - "D[k] = " G ", L[k-1] = " G ", L[k] = " G "\n", - k,gershgorin_tolerance,i, - machine_precision,d,shift,GR, - D[k], (k>0)?L[k-1]:0, (k+1 -#include -#include -#include - -#define min(a,b) ((a)<(b)?(a):(b)) - -void mergesort(double *X, int n, int *index) -{ - for(int i=0;i=r) continue; /* If right list is empty, sequence is already sorted */ - -#if TEST_SORT - printf("Step %d:%d: Merging pre-sorted segments: [", step,i); - for(int k=i;k -#include - -#define FLOAT_TYPE float64 - -#if FLOAT_TYPE==float64 - #define real_t double - #define machine_precision DBL_EPSILON - #define PRINTF_G "%g" - -#elif FLOAT_TYPE==float80 - #define real_t long double - #define PRINTF_G "%Lg" - #define creal(x) creall(x) - #define cimag(x) cimagl(x) - - #define conj conjl - #define sqrt sqrtl - #define fabs fabsl - #define machine_precision LDBL_EPSILON - -#elif FLOAT_TYPE==float32 - #define real_t float - #define PRINTF_G "%g" - #define creal(x) crealf(x) - #define cimag(x) cimagf(x) - - #define conj conjf - #define sqrt sqrtf - #define fabs fabsf - #define machine_precision FLT_EPSILON -#endif - -#ifdef __cplusplus - #include - typedef std::complex complex_t; -#else -typedef real_t complex complex_t; -#endif - - -#define INLINE inline __attribute__((always_inline)) - - -typedef struct { - real_t value[2]; -} real_pair; - -typedef complex_t scalar; From 6d7377a8f310cdb7a0fdae30633329174761d64d Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Tue, 24 Feb 2026 13:51:31 +0100 Subject: [PATCH 04/15] Transport J Avery eigen "lib" to mcstas-comps/share and %include --- mcstas-comps/contrib/Phonon_BvK_PG.comp | 5 +- .../Test_Phonon_BvK_PG.instr | 4 +- mcstas-comps/share/eigen-types.h | 50 ++ mcstas-comps/share/eigen.c | 451 ++++++++++++++++++ mcstas-comps/share/eigen.h | 87 ++++ 5 files changed, 592 insertions(+), 5 deletions(-) create mode 100644 mcstas-comps/share/eigen-types.h create mode 100644 mcstas-comps/share/eigen.c create mode 100644 mcstas-comps/share/eigen.h diff --git a/mcstas-comps/contrib/Phonon_BvK_PG.comp b/mcstas-comps/contrib/Phonon_BvK_PG.comp index bb2812190..dbbd9b0f9 100644 --- a/mcstas-comps/contrib/Phonon_BvK_PG.comp +++ b/mcstas-comps/contrib/Phonon_BvK_PG.comp @@ -75,8 +75,8 @@ OUTPUT PARAMETERS (q_x, q_y, q_z) SHARE %{ -//%include "eigen.h" -//%include "eigen.c" +%include "eigen.h" +%include "eigen.c" #ifndef PHONON_SIMPLE #define PHONON_SIMPLE $Revision$ // KL comment 040125: what is this used for? #pragma acc routine @@ -166,7 +166,6 @@ void bubblesort (int size , double* inputarray , int* index) #include #include #include -#include "eigen.h" extern FILE *matrix_test_file, *parms_test_file; extern int matrix_test_count; diff --git a/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/Test_Phonon_BvK_PG.instr b/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/Test_Phonon_BvK_PG.instr index 187dcef05..aafff80bd 100644 --- a/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/Test_Phonon_BvK_PG.instr +++ b/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/Test_Phonon_BvK_PG.instr @@ -31,8 +31,8 @@ *******************************************************************************/ 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 " +//SHELL " cc -c eigen.c " +//DEPENDENCY " eigen.o " DECLARE %{ diff --git a/mcstas-comps/share/eigen-types.h b/mcstas-comps/share/eigen-types.h new file mode 100644 index 000000000..7a3600e58 --- /dev/null +++ b/mcstas-comps/share/eigen-types.h @@ -0,0 +1,50 @@ +#pragma once +#include +#include + +#define FLOAT_TYPE float64 + +#if FLOAT_TYPE==float64 + #define real_t double + #define machine_precision DBL_EPSILON + #define PRINTF_G "%g" + +#elif FLOAT_TYPE==float80 + #define real_t long double + #define PRINTF_G "%Lg" + #define creal(x) creall(x) + #define cimag(x) cimagl(x) + + #define conj conjl + #define sqrt sqrtl + #define fabs fabsl + #define machine_precision LDBL_EPSILON + +#elif FLOAT_TYPE==float32 + #define real_t float + #define PRINTF_G "%g" + #define creal(x) crealf(x) + #define cimag(x) cimagf(x) + + #define conj conjf + #define sqrt sqrtf + #define fabs fabsf + #define machine_precision FLT_EPSILON +#endif + +#ifdef __cplusplus + #include + typedef std::complex complex_t; +#else +typedef real_t complex complex_t; +#endif + + +#define INLINE inline __attribute__((always_inline)) + + +typedef struct { + real_t value[2]; +} real_pair; + +typedef complex_t scalar; diff --git a/mcstas-comps/share/eigen.c b/mcstas-comps/share/eigen.c new file mode 100644 index 000000000..db990ac9f --- /dev/null +++ b/mcstas-comps/share/eigen.c @@ -0,0 +1,451 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +/* TODO: Tridiagonal input/outpout in (3,m) storage. */ +/* TODO: Symmetric tridiagonal input/outpout in (2,m) storage. */ +/* TODO: Add pivot to handle ill-conditioned matrices */ +//#include "eigen.h" + +#define G PRINTF_G +#define ROWS 0 +#define COLS 1 + +/* TOOD: Generic matrix functions to separate file */ +INLINE real_t vector_norm2(const scalar *x, int n){ + real_t sum = 0; + for(int i=0;i numerical_zero) // Only process if subdiagonal element is not already zero. + { + a[0] = D[k]; a[1] = L[k]; // a = T[k:k+2,k] is the vector of nonzeros in kth subdiagonal column. + + real_t anorm = sqrt(a[0]*a[0] + a[1]*a[1]); + + // // Udrullet + // // reflection_vector(a,anorm,v); + v[0] = D[k]; v[1] = L[k]; + real_t alpha = -copysign(anorm,a[0]); // Koster ingenting + v[0] -= alpha; + + real_t vnorm = sqrt(v[0]*v[0]+v[1]*v[1]); + real_t norm_inv = 1/vnorm; /* Normalize */ + v[0] *= norm_inv; v[1] *= norm_inv; + + Vout[2*k] = v[0]; Vout[2*k+1] = v[1]; + + // // Udrullet + // // apply_reflection(T({k,k+2},{k,k+3}),v); + // // if(k+1=0;k--){ + // We start by targeting the (n,n)-eigenvalue, and gradually + // deflate, working on smaller and smaller submatrices. + real_t d = D[k]; // d = T(k,k) + real_t shift = d; + + // The Gershgorin disk radius is defined by just the row-sums of + // absolute off-diagonal elements, since T is symmetric. As T is + // tridiagonal, only T(k,k-1),T(k,k), and T(k,k+1) are nonzero. + // Thus, the k'th Gershgorin radius is just |T(k,k-1)| + + // |T(k,k+1)| = |T(k,k-1)| + |T(k+1,k)| = |L[k-1]|+|L[k]|. + int i=0; + + real_t GR = (k>0?fabs(L[k-1]):0)+fabs(L[k]); // Initial Gershgorin radius + int not_done = 2; + while(not_done > 0){ + i++; + T_QTQ(k+1, D,L, D,L, V, shift); + if(Q!=0) apply_real_reflections(V,k,Qhat,n); + + GR = (k>0?fabs(L[k-1]):0)+(k+10){ + real_pair l = eigvalsh2x2(D[k-1],L[k-1], /* T[(k-1):k, (k-1):k] */ + L[k-1],D[k] ); + real_t l0 = l.value[0], l1 = l.value[1]; + + shift = fabs(l0-d) < fabs(l1-d)? l0 : l1; // Pick closest eigenvalue + } else + shift = D[k]; + + if(GR <= gershgorin_tolerance) not_done--; + if(i>max_iterations && GR > gershgorin_tolerance){ + fprintf(stderr,"Cannot converge eigenvalue %d to tolerance " G + " in %d iterations using machine precision " G " (d=" G ", shift=" G ", GR=" G ")\n" + "D[k] = " G ", L[k-1] = " G ", L[k] = " G "\n", + k,gershgorin_tolerance,i, + machine_precision,d,shift,GR, + D[k], (k>0)?L[k-1]:0, (k+1 +#include + +#define FLOAT_TYPE float64 + +#if FLOAT_TYPE==float64 + #define real_t double + #define machine_precision DBL_EPSILON + #define PRINTF_G "%g" + +#elif FLOAT_TYPE==float80 + #define real_t long double + #define PRINTF_G "%Lg" + #define creal(x) creall(x) + #define cimag(x) cimagl(x) + + #define conj conjl + #define sqrt sqrtl + #define fabs fabsl + #define machine_precision LDBL_EPSILON + +#elif FLOAT_TYPE==float32 + #define real_t float + #define PRINTF_G "%g" + #define creal(x) crealf(x) + #define cimag(x) cimagf(x) + + #define conj conjf + #define sqrt sqrtf + #define fabs fabsf + #define machine_precision FLT_EPSILON +#endif + +#ifdef __cplusplus + #include + typedef std::complex complex_t; +#else +typedef real_t complex complex_t; +#endif + + +#define INLINE inline __attribute__((always_inline)) + + +typedef struct { + real_t value[2]; +} real_pair; + +typedef complex_t scalar; + + +void print_vector(const char *name, scalar *a, int l); +void print_matrix(const char *name, scalar *A, int m, int n); + +real_t vector_norm(const scalar *x, int n); +void extract_region(scalar *S, int N, + int i0, int j0, int m, int n, + scalar *D); +real_t max_norm(complex_t *A, int m, int n); +void identity_matrix(scalar *Q, int n); +void real_identity_matrix(real_t *Q, int n); +void matrix_inplace_multiply(scalar *A, const scalar *B, + int m, int n, int q); + +void reflection_vector(/*in*/const scalar *a, const real_t anorm, + /*out*/scalar *v, scalar *sigma, int n); +void apply_reflection(/*in/out*/scalar *A, const scalar *v, + int m, int n, scalar sigma, int transpose); + +void reflect_region(/*in/out*/scalar *A, int N, + int i0, int j0, int m, int n, + const scalar *v, scalar sigma, int cols); + +void apply_real_reflections(const real_t *V, const int n, real_t *Q, const int m); + + +void QHQ(/*in/out*/scalar *A, int n, scalar *Q); +void T_QTQ(const int n, + const real_t *Din, const real_t *Lin, + real_t *Dout, real_t *Lout, real_t *Vout, + real_t shift); + +real_pair eigvalsh2x2(real_t a, real_t b, real_t c, real_t d); + +real_pair eigensystem_hermitian(const scalar *A, const int n, + real_t *lambdas, scalar *Q); + From dec85d45ccc9d70d304fd856bee1cf968881095a Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Tue, 24 Feb 2026 13:53:40 +0100 Subject: [PATCH 05/15] Remove eigen-types.h, was added top of eigen.h --- mcstas-comps/share/eigen-types.h | 50 -------------------------------- 1 file changed, 50 deletions(-) delete mode 100644 mcstas-comps/share/eigen-types.h diff --git a/mcstas-comps/share/eigen-types.h b/mcstas-comps/share/eigen-types.h deleted file mode 100644 index 7a3600e58..000000000 --- a/mcstas-comps/share/eigen-types.h +++ /dev/null @@ -1,50 +0,0 @@ -#pragma once -#include -#include - -#define FLOAT_TYPE float64 - -#if FLOAT_TYPE==float64 - #define real_t double - #define machine_precision DBL_EPSILON - #define PRINTF_G "%g" - -#elif FLOAT_TYPE==float80 - #define real_t long double - #define PRINTF_G "%Lg" - #define creal(x) creall(x) - #define cimag(x) cimagl(x) - - #define conj conjl - #define sqrt sqrtl - #define fabs fabsl - #define machine_precision LDBL_EPSILON - -#elif FLOAT_TYPE==float32 - #define real_t float - #define PRINTF_G "%g" - #define creal(x) crealf(x) - #define cimag(x) cimagf(x) - - #define conj conjf - #define sqrt sqrtf - #define fabs fabsf - #define machine_precision FLT_EPSILON -#endif - -#ifdef __cplusplus - #include - typedef std::complex complex_t; -#else -typedef real_t complex complex_t; -#endif - - -#define INLINE inline __attribute__((always_inline)) - - -typedef struct { - real_t value[2]; -} real_pair; - -typedef complex_t scalar; From 7b0914cfcc098dca79c4087f14ef7cc0502ff021 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Tue, 24 Feb 2026 14:04:38 +0100 Subject: [PATCH 06/15] Remove OUTPUT and DEFINITION pars (McStas 2.x style) --- mcstas-comps/contrib/Phonon_BvK_PG.comp | 2 -- 1 file changed, 2 deletions(-) diff --git a/mcstas-comps/contrib/Phonon_BvK_PG.comp b/mcstas-comps/contrib/Phonon_BvK_PG.comp index dbbd9b0f9..ebd3ef344 100644 --- a/mcstas-comps/contrib/Phonon_BvK_PG.comp +++ b/mcstas-comps/contrib/Phonon_BvK_PG.comp @@ -67,10 +67,8 @@ ******************************************************************************/ DEFINE COMPONENT Phonon_BvK_PG -DEFINITION PARAMETERS () 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) -OUTPUT PARAMETERS (q_x, q_y, q_z) /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ SHARE From 470e7d926401eb35608cf9390d599db473a1f31d Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Tue, 24 Feb 2026 14:14:58 +0100 Subject: [PATCH 07/15] Use mccode-style single-line %include import --- mcstas-comps/contrib/Phonon_BvK_PG.comp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/mcstas-comps/contrib/Phonon_BvK_PG.comp b/mcstas-comps/contrib/Phonon_BvK_PG.comp index ebd3ef344..194b4933d 100644 --- a/mcstas-comps/contrib/Phonon_BvK_PG.comp +++ b/mcstas-comps/contrib/Phonon_BvK_PG.comp @@ -73,8 +73,9 @@ focus_r=0,focus_xw=0,focus_yh=0,focus_aw=0,focus_ah=0,target_x=0, target_y=0, ta SHARE %{ -%include "eigen.h" -%include "eigen.c" + // James Avery eigenvalue solver + %include "eigen" + #ifndef PHONON_SIMPLE #define PHONON_SIMPLE $Revision$ // KL comment 040125: what is this used for? #pragma acc routine From 308ace44ea966697ad44223bb2a784a63684a564 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Tue, 24 Feb 2026 14:15:20 +0100 Subject: [PATCH 08/15] Protecti for multiple imports via defines --- mcstas-comps/share/eigen.c | 5 ++++- mcstas-comps/share/eigen.h | 5 +++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/mcstas-comps/share/eigen.c b/mcstas-comps/share/eigen.c index db990ac9f..34925f1f3 100644 --- a/mcstas-comps/share/eigen.c +++ b/mcstas-comps/share/eigen.c @@ -1,3 +1,6 @@ +#ifndef EIGEN_LIB +#define EIGEN_LIB + #include #include #include @@ -448,4 +451,4 @@ real_pair eigensystem_hermitian(const scalar *A, } - +#endif diff --git a/mcstas-comps/share/eigen.h b/mcstas-comps/share/eigen.h index fdcd5fa3d..a8500d114 100644 --- a/mcstas-comps/share/eigen.h +++ b/mcstas-comps/share/eigen.h @@ -1,3 +1,6 @@ +#ifndef EIGEN_LIB_H +#define EIGEN_LIB_H + #include #include @@ -85,3 +88,5 @@ real_pair eigvalsh2x2(real_t a, real_t b, real_t c, real_t d); real_pair eigensystem_hermitian(const scalar *A, const int n, real_t *lambdas, scalar *Q); + +#endif From bacaad9bab96c2caf5dcc657a62548e26b6d922b Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Tue, 24 Feb 2026 14:24:24 +0100 Subject: [PATCH 09/15] Add headers mentioning James Avery contribution --- mcstas-comps/share/eigen.c | 10 +++++++++- mcstas-comps/share/eigen.h | 12 +++++++++++- 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/mcstas-comps/share/eigen.c b/mcstas-comps/share/eigen.c index 34925f1f3..ab49f05c3 100644 --- a/mcstas-comps/share/eigen.c +++ b/mcstas-comps/share/eigen.c @@ -1,3 +1,11 @@ +/************************************************************ + * McStas Eigensolver library * + * Contributed by James Emil Avery, 202x * + * Department of Computer Science, University of Copenhagen * + ************************************************************/ + +/* Original file "eigen.c": */ + #ifndef EIGEN_LIB #define EIGEN_LIB @@ -13,7 +21,7 @@ /* TODO: Tridiagonal input/outpout in (3,m) storage. */ /* TODO: Symmetric tridiagonal input/outpout in (2,m) storage. */ /* TODO: Add pivot to handle ill-conditioned matrices */ -//#include "eigen.h" +// #include "eigen.h" #define G PRINTF_G #define ROWS 0 diff --git a/mcstas-comps/share/eigen.h b/mcstas-comps/share/eigen.h index a8500d114..47e6a3b0a 100644 --- a/mcstas-comps/share/eigen.h +++ b/mcstas-comps/share/eigen.h @@ -1,6 +1,13 @@ +/************************************************************ + * McStas Eigensolver library * + * Contributed by James Emil Avery, 202x * + * Department of Computer Science, University of Copenhagen * + ************************************************************/ + #ifndef EIGEN_LIB_H #define EIGEN_LIB_H +/* Original file "types.h": */ #include #include @@ -51,6 +58,9 @@ typedef struct { typedef complex_t scalar; +/* End of "types.h" */ +/* ------------------------ */ +/* Original file "eigen.h": */ void print_vector(const char *name, scalar *a, int l); void print_matrix(const char *name, scalar *A, int m, int n); @@ -88,5 +98,5 @@ real_pair eigvalsh2x2(real_t a, real_t b, real_t c, real_t d); real_pair eigensystem_hermitian(const scalar *A, const int n, real_t *lambdas, scalar *Q); - +/* End of eigen.h */ #endif From 4eb771facacef4c04a7a01d44079a23aea8bb2c2 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Tue, 24 Feb 2026 14:33:21 +0100 Subject: [PATCH 10/15] Mark component NOACC (execute CPU only for now) --- mcstas-comps/contrib/Phonon_BvK_PG.comp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/mcstas-comps/contrib/Phonon_BvK_PG.comp b/mcstas-comps/contrib/Phonon_BvK_PG.comp index 194b4933d..e81ad6ad4 100644 --- a/mcstas-comps/contrib/Phonon_BvK_PG.comp +++ b/mcstas-comps/contrib/Phonon_BvK_PG.comp @@ -69,7 +69,8 @@ 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) -/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ +NOACC +// The component is currently "NOACC" only. SHARE %{ From 6351866c0f22ee2d039544c817dbba1dfb1d2929 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 25 Feb 2026 14:32:15 +0100 Subject: [PATCH 11/15] Use mccode-complex-lib --- mcstas-comps/contrib/Phonon_BvK_PG.comp | 1 + 1 file changed, 1 insertion(+) diff --git a/mcstas-comps/contrib/Phonon_BvK_PG.comp b/mcstas-comps/contrib/Phonon_BvK_PG.comp index e81ad6ad4..4a9d80cb4 100644 --- a/mcstas-comps/contrib/Phonon_BvK_PG.comp +++ b/mcstas-comps/contrib/Phonon_BvK_PG.comp @@ -74,6 +74,7 @@ NOACC SHARE %{ + %include "mccode-complex-lib" // James Avery eigenvalue solver %include "eigen" From 578de7c8e360e7eb170f2080c07c737c20ee948e Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 25 Feb 2026 15:24:06 +0100 Subject: [PATCH 12/15] Let J Avery solver use mccode-complex-lib cdouble type --- mcstas-comps/share/eigen.c | 126 ++++++++++++++++++++----------------- mcstas-comps/share/eigen.h | 83 ++++++++---------------- 2 files changed, 94 insertions(+), 115 deletions(-) diff --git a/mcstas-comps/share/eigen.c b/mcstas-comps/share/eigen.c index ab49f05c3..67957f47c 100644 --- a/mcstas-comps/share/eigen.c +++ b/mcstas-comps/share/eigen.c @@ -28,29 +28,29 @@ #define COLS 1 /* TOOD: Generic matrix functions to separate file */ -INLINE real_t vector_norm2(const scalar *x, int n){ +INLINE real_t vector_norm2(const cdouble*x, int n){ real_t sum = 0; - for(int i=0;i -#include -#define FLOAT_TYPE float64 - -#if FLOAT_TYPE==float64 - #define real_t double - #define machine_precision DBL_EPSILON - #define PRINTF_G "%g" - -#elif FLOAT_TYPE==float80 - #define real_t long double - #define PRINTF_G "%Lg" - #define creal(x) creall(x) - #define cimag(x) cimagl(x) - - #define conj conjl - #define sqrt sqrtl - #define fabs fabsl - #define machine_precision LDBL_EPSILON - -#elif FLOAT_TYPE==float32 - #define real_t float +// PW: Non-float64 logic suppressed. + #define real_t double + #define scalar cdouble + #define machine_precision DBL_EPSILON #define PRINTF_G "%g" - #define creal(x) crealf(x) - #define cimag(x) cimagf(x) - - #define conj conjf - #define sqrt sqrtf - #define fabs fabsf - #define machine_precision FLT_EPSILON -#endif -#ifdef __cplusplus - #include - typedef std::complex complex_t; +#ifndef _MSC_EXTENSIONS +#define INLINE inline __attribute__((always_inline)) #else -typedef real_t complex complex_t; +#define INLINE inline #endif - -#define INLINE inline __attribute__((always_inline)) - - typedef struct { real_t value[2]; } real_pair; -typedef complex_t scalar; /* End of "types.h" */ /* ------------------------ */ /* Original file "eigen.h": */ -void print_vector(const char *name, scalar *a, int l); -void print_matrix(const char *name, scalar *A, int m, int n); +void print_vector(const char *name, cdouble *a, int l); +void print_matrix(const char *name, cdouble *A, int m, int n); -real_t vector_norm(const scalar *x, int n); -void extract_region(scalar *S, int N, +real_t vector_norm(const cdouble *x, int n); +void extract_region(cdouble *S, int N, int i0, int j0, int m, int n, - scalar *D); -real_t max_norm(complex_t *A, int m, int n); -void identity_matrix(scalar *Q, int n); -void real_identity_matrix(real_t *Q, int n); -void matrix_inplace_multiply(scalar *A, const scalar *B, + cdouble *D); +double max_norm(cdouble* A, int m, int n); +void identity_matrix(cdouble *Q, int n); +void real_identity_matrix(double *Q, int n); +void matrix_inplace_multiply(cdouble *A, cdouble *B, int m, int n, int q); -void reflection_vector(/*in*/const scalar *a, const real_t anorm, - /*out*/scalar *v, scalar *sigma, int n); -void apply_reflection(/*in/out*/scalar *A, const scalar *v, - int m, int n, scalar sigma, int transpose); +void reflection_vector(/*in*/cdouble *a, double anorm, + /*out*/cdouble *v, cdouble *sigma, int n); +void apply_reflection(/*in/out*/cdouble *A, const cdouble *v, + int m, int n, cdouble sigma, int transpose); -void reflect_region(/*in/out*/scalar *A, int N, +void reflect_region(/*in/out*/cdouble *A, int N, int i0, int j0, int m, int n, - const scalar *v, scalar sigma, int cols); + const cdouble *v, cdouble sigma, int cols); -void apply_real_reflections(const real_t *V, const int n, real_t *Q, const int m); +void apply_real_reflections(double *V, int n, double *Q, int m); -void QHQ(/*in/out*/scalar *A, int n, scalar *Q); +void QHQ(/*in/out*/cdouble *A, int n, cdouble *Q); void T_QTQ(const int n, const real_t *Din, const real_t *Lin, real_t *Dout, real_t *Lout, real_t *Vout, @@ -95,8 +64,8 @@ void T_QTQ(const int n, real_pair eigvalsh2x2(real_t a, real_t b, real_t c, real_t d); -real_pair eigensystem_hermitian(const scalar *A, const int n, - real_t *lambdas, scalar *Q); +real_pair eigensystem_hermitian(const cdouble *A, int n, + real_t *lambdas, cdouble *Q); /* End of eigen.h */ #endif From 0b98c9e124510cde3846840a282144071eb33d4e Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 25 Feb 2026 21:01:31 +0100 Subject: [PATCH 13/15] Apply clangformat for reability/structure --- mcstas-comps/contrib/Phonon_BvK_PG.comp | 2031 +++++++++++------------ 1 file changed, 1008 insertions(+), 1023 deletions(-) diff --git a/mcstas-comps/contrib/Phonon_BvK_PG.comp b/mcstas-comps/contrib/Phonon_BvK_PG.comp index 4a9d80cb4..45c896122 100644 --- a/mcstas-comps/contrib/Phonon_BvK_PG.comp +++ b/mcstas-comps/contrib/Phonon_BvK_PG.comp @@ -78,863 +78,859 @@ SHARE // James Avery eigenvalue solver %include "eigen" -#ifndef PHONON_SIMPLE -#define PHONON_SIMPLE $Revision$ // KL comment 040125: what is this used for? -#pragma acc routine + #ifndef PHONON_SIMPLE + #define PHONON_SIMPLE $Revision$ // KL comment 040125: what is this used for? + #pragma acc routine -// ---------------- THINGS THAT COULD BE GENERAL FOR McStas -------------- + // ---------------- THINGS THAT COULD BE GENERAL FOR McStas -------------- -#define T2E (1/11.605) // Kelvin to meV converter -#define pi 3.14159265358979323846264338327950288419716939937510L -#define sqrt3 1.73205080756887729352744634150587236694280525381038L -#define THz2meV 4.13566 // THz to meV converter + #define T2E (1/11.605) // Kelvin to meV converter + #define pi 3.14159265358979323846264338327950288419716939937510L + #define sqrt3 1.73205080756887729352744634150587236694280525381038L + #define THz2meV 4.13566 // THz to meV converter -#define Da2kg 1.66054e-27 //Dalton to kg converter -#define dyn2N 1e-3 // dyn/cm to N/m converter -#define hbar 1.0546E-34 -#define mn (1.0087*Da2kg) + #define Da2kg 1.66054e-27 //Dalton to kg converter + #define dyn2N 1e-3 // dyn/cm to N/m converter + #define hbar 1.0546E-34 + #define mn (1.0087*Da2kg) -double nbose(double omega, double T) /* Give it another name ?? */ + double + nbose (double omega, double T) /* Give it another name ?? */ { double nb; - nb= (omega>0) ? 1+1/(exp(omega/(T*T2E))-1) : 1/(exp(-omega/(T*T2E))-1); + nb = (omega > 0) ? 1 + 1 / (exp (omega / (T * T2E)) - 1) : 1 / (exp (-omega / (T * T2E)) - 1); return nb; } -#undef T2E + #undef T2E -/* Routine types from Numerical Recipies book */ -#define UNUSED (-1.11e30) -#define MAXRIDD 60 + /* Routine types from Numerical Recipies book */ + #define UNUSED (-1.11e30) + #define MAXRIDD 60 -void fatalerror_cpu(char *s) - { - fprintf(stderr,"%s \n",s); - exit(1); + void + fatalerror_cpu (char* s) { + fprintf (stderr, "%s \n", s); + exit (1); } - void nrerror(const char *error_text){ - printf("Numerical Recipes run-time error ...\n"); - printf("%s\n",error_text); - printf("... now exiting to system.\n"); - exit(1); -} + void + nrerror (const char* error_text) { + printf ("Numerical Recipes run-time error ...\n"); + printf ("%s\n", error_text); + printf ("... now exiting to system.\n"); + exit (1); + } -#pragma acc routine -void fatalerror(char *s) - { - #ifndef OPENACC - fatalerror_cpu(s); - #endif + #pragma acc routine + void + fatalerror (char* s) { + #ifndef OPENACC + fatalerror_cpu (s); + #endif } #pragma acc routine -void bubblesort (int size , double* inputarray , int* index) -//this function modifies the inputarray by bubble sorting, and also modifies the index array as the -//index after the sorting. The index array should be initialized as index[] = {0,1,2,3,4,5,6,7...} -{ + void + bubblesort (int size, double* inputarray, int* index) + // this function modifies the inputarray by bubble sorting, and also modifies the index array as the + // index after the sorting. The index array should be initialized as index[] = {0,1,2,3,4,5,6,7...} + { int tempindex = 0; int i; int j; - for (i = 0 ; i < size-1 ; i++) - { - for (j = 0 ; j < size-1 ; j++) - { - if (inputarray[index[j]] > inputarray[index[j+1]]) - { - tempindex = index[j]; - index[j] = index[j+1]; - index[j+1] = tempindex; - } + for (i = 0; i < size - 1; i++) { + for (j = 0; j < size - 1; j++) { + if (inputarray[index[j]] > inputarray[index[j + 1]]) { + tempindex = index[j]; + index[j] = index[j + 1]; + index[j + 1] = tempindex; } - } + } + } return; -} - -// ---------------- END GENERAL FOR McStas ---------------------------------- - -// ---------------- THINGS THAT ARE SPECIFIC FOR THE PG COMPONENT ------------ - -#define SMALL_NUMBER 1E-6 -#define V_HIGH 8000 // Highest velocity to search, corresponding to 320 meV, far above the to of the PG phonon band (202 meV along main axes) -#define MATRIX_TEST 0 // Change to 1 to write matrices to disk -#define TIME_EIGENSYSTEMS 1 - -#define DIM 12 // dimensionality of the phonon eigenvectors (4 atoms in the unit cell times 3 dimensions) -#define DV 0.001 // Velocity change used for numerical derivative [m/s] CHECK if this should be smaller - -#include -#include -#include -#include - -extern FILE *matrix_test_file, *parms_test_file; -extern int matrix_test_count; -int mode, verbose, shape; - -#define a_latt 2.461 //PG lattice constant in AA; Nicklow1972 gives 2.45 (WHICH PAGE?) -#define c_latt 6.708 // PG lattice constant in AA; Nicklow1972 gives 6.70 -#define b_length 6.646 // PG scattering legth in fm - -double M = 12.011 * Da2kg; //atomic mass of C - -#define K_l1 3.62e+5 // Force constant longitudinal nn1 - in (a,b) plane; here units of dyn -#define K_t1 1.99e+5 // Force constant transverse nn1 - in (a,b) plane -#define K_l2 1.33e+5 // Force constant longitudinal nn2 - in (a,b) plane -#define K_t2 -0.520e+5 // Force constant transverse nn2 - in (a,b) plane -#define K_l3 -0.037e+5 // Force constant longitudinal nn3 - in (a,b) plane -#define K_t3 0.288e+5 // Force constant transverse nn3 - in (a,b) plane -#define K_l4 0.058e+5 // Force constant longitudinal nn4 - along c -#define K_t4 0.0077e+5 // Force constant transverse nn4 - along c - - //Lattice basis - - // PG is hexagonal close packed with 4 atoms per cell - // Coordinates from Nicklow1972, Fig. 7: Atoms A, B, C, D - double Delta[4*3] = {0 , 0 , 0 , a_latt/(2*sqrt3) , a_latt/2 , 0 , 0 , 0 , c_latt/2 , -a_latt/(2*sqrt3) , a_latt/2 , c_latt/2}; + } - // Lattice vectors from paper fig. 7 - double avec[3] = {a_latt*sqrt3/2 , a_latt*0.5 , 0}; // THIS IS NEVER USED ??!! - double bvec[3] = {a_latt*(-sqrt3/2) , a_latt*0.5 , 0}; - double cvec[3] = {0 , 0 , c_latt}; + // ---------------- END GENERAL FOR McStas ---------------------------------- - // reciprocal lattice vectors - double astar = 4*PI/(sqrt3*a_latt); // length of reciprocal lattice vector a* - double cstar = 2*PI/c_latt; // length of reciprocal lattice vector c* + // ---------------- THINGS THAT ARE SPECIFIC FOR THE PG COMPONENT ------------ - // Rotation matrices - // m(12*i+j) = M(i,j) - double Rot120[3][3] = {{-0.5 , sqrt3/2 , 0} , {-sqrt3/2 , -0.5 , 0} , {0 , 0 , 1}}; - double Rot120_flat[9] = {-0.5 , sqrt3/2 , 0 , -sqrt3/2 , -0.5 , 0 , 0 , 0 , 1}; - double Rot60[3][3] = {{0.5 , sqrt3/2 , 0} , {-sqrt3/2 , 0.5 , 0} , {0 , 0 , 1}}; - double Rot120rev[3][3] = {{-0.5 , -sqrt3/2 , 0} , {sqrt3/2 , -0.5 , 0} , {0 , 0 , 1}}; //rotate -120 - double Rot60rev[3][3] = {{0.5 , -sqrt3/2 , 0} , {sqrt3/2 , 0.5 , 0} , {0 , 0 , 1}}; //rotate -60 - double Rot180[3][3] = {{-1, 0, 0}, {0, -1, 0}, {0, 0, 1}}; // rotate 180 or -180 -// double Rot5[3][3] = {{0.5 , sqrt3/2 , 0} , {-sqrt3/2 , 0.5 , 0} , {0 , 0 , 1}}; //rotate 60 -// double Rot5rev[3][3] = {{0.5 , -sqrt3/2 , 0} , {sqrt3/2 , 0.5 , 0} , {0 , 0 , 1}}; //rotate -60 + #define SMALL_NUMBER 1E-6 + #define V_HIGH 8000 // Highest velocity to search, corresponding to 320 meV, far above the to of the PG phonon band (202 meV along main axes) + #define MATRIX_TEST 0 // Change to 1 to write matrices to disk + #define TIME_EIGENSYSTEMS 1 - //1st neighbour + #define DIM 12 // dimensionality of the phonon eigenvectors (4 atoms in the unit cell times 3 dimensions) + #define DV 0.001 // Velocity change used for numerical derivative [m/s] CHECK if this should be smaller -// double r_j1[3][3] = {{-a_latt/sqrt3 , 0 , 0} , {a_latt*0.5/sqrt3 , -a_latt*0.5 , 0} , {a_latt*0.5/sqrt3 , a_latt*0.5 , 0}}; // This holds from atoms A and D -// double r_j2[3][3] = {{a_latt/sqrt3 , 0 , 0} , {-a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {-a_latt*0.5/sqrt3 , -a_latt*0.5 , 0}}; // This holds from atoms B and C + #include + #include + #include + #include - double Phi_nn1[3][3] = {{K_l1*dyn2N , 0 , 0} , {0 , K_t1*dyn2N , 0} , {0 , 0 , K_t1*dyn2N}}; - double Phi1[3][3] = {{K_l1*dyn2N , 0 , 0} , {0 , K_t1*dyn2N , 0} , {0 , 0 , K_t1*dyn2N}}; //Phi1 = Phi_nn1 - double Phi2[3][3]; - double Phi3[3][3]; + extern FILE *matrix_test_file, *parms_test_file; + extern int matrix_test_count; + int mode, verbose, shape; - //2nd neighbour -// double r_j11[9][3] = {{a_latt*(-1/sqrt3) , 0 , 0} , {a_latt*0.5/sqrt3 , a_latt*(-0.5) , 0} , {a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {0 , a_latt , 0} , {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , 0}};//r_j1 + r_add -// double r_j21[9][3] = {{-a_latt*(-1/sqrt3) , -a_latt*0 , 0} , {-a_latt*0.5/sqrt3 , -a_latt*(-0.5) , 0} , {-a_latt*0.5/sqrt3 , -a_latt*0.5 , 0} , {0 , a_latt , 0} , {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , 0}};//r_j2 + r_add + #define a_latt 2.461 //PG lattice constant in AA; Nicklow1972 gives 2.45 (WHICH PAGE?) + #define c_latt 6.708 // PG lattice constant in AA; Nicklow1972 gives 6.70 + #define b_length 6.646 // PG scattering legth in fm - double Phi_nn2[3][3] = {{K_t2*dyn2N , 0 , 0} , {0 , K_l2*dyn2N , 0} , {0 , 0 , K_t2*dyn2N}}; - double Phi4[3][3] = {{K_t2*dyn2N , 0 , 0} , {0 , K_l2*dyn2N , 0} , {0 , 0 , K_t2*dyn2N}}; - double Phi5[3][3]; - double Phi6[3][3]; - double Phi7[3][3]; - double Phi8[3][3]; - double Phi9[3][3]; + double M = 12.011 * Da2kg; // atomic mass of C - //3rd neighbour -// double r_j12[12][3] = {{a_latt*(-1/sqrt3) , 0 , 0} , {a_latt*0.5/sqrt3 , a_latt*(-0.5) , 0} , {a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {0 , a_latt , 0} , {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , 0} , {2*a_latt/sqrt3 , 0 , 0} , {-a_latt/sqrt3 , a_latt , 0} , {-a_latt/sqrt3 , -a_latt , 0}};//r_j11 + r_add2 - double r_j2[12][3] = {{a_latt/sqrt3 , 0 , 0} , {-a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {-a_latt*0.5/sqrt3 , -a_latt*0.5 , 0} , {0 , a_latt , 0} , {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , 0} , {-2*a_latt/sqrt3 , 0 , 0} , {a_latt/sqrt3 , -a_latt , 0} , {a_latt/sqrt3 , a_latt , 0}};//r_j21 - r_add2 + #define K_l1 3.62e+5 // Force constant longitudinal nn1 - in (a,b) plane; here units of dyn + #define K_t1 1.99e+5 // Force constant transverse nn1 - in (a,b) plane + #define K_l2 1.33e+5 // Force constant longitudinal nn2 - in (a,b) plane + #define K_t2 -0.520e+5 // Force constant transverse nn2 - in (a,b) plane + #define K_l3 -0.037e+5 // Force constant longitudinal nn3 - in (a,b) plane + #define K_t3 0.288e+5 // Force constant transverse nn3 - in (a,b) plane + #define K_l4 0.058e+5 // Force constant longitudinal nn4 - along c + #define K_t4 0.0077e+5 // Force constant transverse nn4 - along c - double Phi_nn3[3][3] = {{K_l3*dyn2N , 0 , 0} , {0 , K_t3*dyn2N , 0} , {0 , 0 , K_t3*dyn2N}}; + // Lattice basis - // c2 = 2/sqrt(6); - // c1 = 1/sqrt(6); - // c0 = 1/sqrt(2); - // c3 = 1/sqrt(3); + // PG is hexagonal close packed with 4 atoms per cell + // Coordinates from Nicklow1972, Fig. 7: Atoms A, B, C, D + double Delta[4 * 3] = { 0, 0, 0, a_latt / (2 * sqrt3), a_latt / 2, 0, 0, 0, c_latt / 2, -a_latt / (2 * sqrt3), a_latt / 2, c_latt / 2 }; - double Phi10[3][3] = {{K_l3*dyn2N , 0 , 0} , {0 , K_t3*dyn2N , 0} , {0 , 0 , K_t3*dyn2N}}; //Phi_nn3 - double Phi11[3][3]; - double Phi12[3][3]; + // Lattice vectors from paper fig. 7 + double avec[3] = { a_latt * sqrt3 / 2, a_latt * 0.5, 0 }; // THIS IS NEVER USED ??!! + double bvec[3] = { a_latt * (-sqrt3 / 2), a_latt * 0.5, 0 }; + double cvec[3] = { 0, 0, c_latt }; - //4th neighbour - double r_j13[14][3] = {{-a_latt/sqrt3 , 0 , 0} , {a_latt*0.5/sqrt3 , -a_latt*0.5 , 0} , {a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {0 , a_latt , 0} , {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , 0} , {2*a_latt/sqrt3 , 0 , 0} , {-a_latt/sqrt3 , a_latt , 0} , {-a_latt/sqrt3 , -a_latt , 0} , {0 , 0 , c_latt/2} , {0 , 0 , -c_latt/2}};//r_j12 + r_add3 Sublattice B and D do not have this coupling + // reciprocal lattice vectors + double astar = 4 * PI / (sqrt3 * a_latt); // length of reciprocal lattice vector a* + double cstar = 2 * PI / c_latt; // length of reciprocal lattice vector c* - double Phi_nn4[3][3] = {{K_t4*dyn2N , 0 , 0} , {0 , K_t4*dyn2N , 0} , {0 , 0 , K_l4*dyn2N}}; - - double Phi13[3][3] = {{K_t4*dyn2N , 0 , 0} , {0 , K_t4*dyn2N , 0} , {0 , 0 , K_l4*dyn2N}}; //Phi_nn4; - double Phi14[3][3] = {{K_t4*dyn2N , 0 , 0} , {0 , K_t4*dyn2N , 0} , {0 , 0 , K_l4*dyn2N}}; //Phi_nn4; - -// ----------------- END SPECIFIC FOR THE PG COPONENT -------------------------- + // Rotation matrices + // m(12*i+j) = M(i,j) + double Rot120[3][3] = { { -0.5, sqrt3 / 2, 0 }, { -sqrt3 / 2, -0.5, 0 }, { 0, 0, 1 } }; + double Rot120_flat[9] = { -0.5, sqrt3 / 2, 0, -sqrt3 / 2, -0.5, 0, 0, 0, 1 }; + double Rot60[3][3] = { { 0.5, sqrt3 / 2, 0 }, { -sqrt3 / 2, 0.5, 0 }, { 0, 0, 1 } }; + double Rot120rev[3][3] = { { -0.5, -sqrt3 / 2, 0 }, { sqrt3 / 2, -0.5, 0 }, { 0, 0, 1 } }; // rotate -120 + double Rot60rev[3][3] = { { 0.5, -sqrt3 / 2, 0 }, { sqrt3 / 2, 0.5, 0 }, { 0, 0, 1 } }; // rotate -60 + double Rot180[3][3] = { { -1, 0, 0 }, { 0, -1, 0 }, { 0, 0, 1 } }; // rotate 180 or -180 + // double Rot5[3][3] = {{0.5 , sqrt3/2 , 0} , {-sqrt3/2 , 0.5 , 0} , {0 , 0 , 1}}; //rotate 60 + // double Rot5rev[3][3] = {{0.5 , -sqrt3/2 , 0} , {sqrt3/2 , 0.5 , 0} , {0 , 0 , 1}}; //rotate -60 + + // 1st neighbour + + // double r_j1[3][3] = {{-a_latt/sqrt3 , 0 , 0} , {a_latt*0.5/sqrt3 , -a_latt*0.5 , 0} , {a_latt*0.5/sqrt3 , a_latt*0.5 , 0}}; // This holds from atoms A and + // D double r_j2[3][3] = {{a_latt/sqrt3 , 0 , 0} , {-a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {-a_latt*0.5/sqrt3 , -a_latt*0.5 , 0}}; // This holds from atoms B + // and C + + double Phi_nn1[3][3] = { { K_l1 * dyn2N, 0, 0 }, { 0, K_t1* dyn2N, 0 }, { 0, 0, K_t1* dyn2N } }; + double Phi1[3][3] = { { K_l1 * dyn2N, 0, 0 }, { 0, K_t1* dyn2N, 0 }, { 0, 0, K_t1* dyn2N } }; // Phi1 = Phi_nn1 + double Phi2[3][3]; + double Phi3[3][3]; + + // 2nd neighbour + // double r_j11[9][3] = {{a_latt*(-1/sqrt3) , 0 , 0} , {a_latt*0.5/sqrt3 , a_latt*(-0.5) , 0} , {a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {0 , a_latt , 0} , + // {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , + // 0}};//r_j1 + r_add double r_j21[9][3] = {{-a_latt*(-1/sqrt3) , -a_latt*0 , 0} , {-a_latt*0.5/sqrt3 , -a_latt*(-0.5) , 0} , {-a_latt*0.5/sqrt3 , + // -a_latt*0.5 , 0} , {0 , a_latt , 0} , {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , + // -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , 0}};//r_j2 + r_add + + double Phi_nn2[3][3] = { { K_t2 * dyn2N, 0, 0 }, { 0, K_l2* dyn2N, 0 }, { 0, 0, K_t2* dyn2N } }; + double Phi4[3][3] = { { K_t2 * dyn2N, 0, 0 }, { 0, K_l2* dyn2N, 0 }, { 0, 0, K_t2* dyn2N } }; + double Phi5[3][3]; + double Phi6[3][3]; + double Phi7[3][3]; + double Phi8[3][3]; + double Phi9[3][3]; + + // 3rd neighbour + // double r_j12[12][3] = {{a_latt*(-1/sqrt3) , 0 , 0} , {a_latt*0.5/sqrt3 , a_latt*(-0.5) , 0} , {a_latt*0.5/sqrt3 , a_latt*0.5 , 0} , {0 , a_latt , 0} , + // {-a_latt*sqrt3/2 , a_latt/2 , 0} , {-a_latt*sqrt3/2 , -a_latt/2 , 0} , {0 , -a_latt , 0} , {a_latt*sqrt3/2 , -a_latt/2 , 0} , {a_latt*sqrt3/2 , a_latt/2 , + // 0} , {2*a_latt/sqrt3 , 0 , 0} , {-a_latt/sqrt3 , a_latt , 0} , {-a_latt/sqrt3 , -a_latt , 0}};//r_j11 + r_add2 + double r_j2[12][3] = { { a_latt / sqrt3, 0, 0 }, + { -a_latt * 0.5 / sqrt3, a_latt * 0.5, 0 }, + { -a_latt * 0.5 / sqrt3, -a_latt * 0.5, 0 }, + { 0, a_latt, 0 }, + { -a_latt * sqrt3 / 2, a_latt / 2, 0 }, + { -a_latt * sqrt3 / 2, -a_latt / 2, 0 }, + { 0, -a_latt, 0 }, + { a_latt * sqrt3 / 2, -a_latt / 2, 0 }, + { a_latt * sqrt3 / 2, a_latt / 2, 0 }, + { -2 * a_latt / sqrt3, 0, 0 }, + { a_latt / sqrt3, -a_latt, 0 }, + { a_latt / sqrt3, a_latt, 0 } }; // r_j21 - r_add2 + + double Phi_nn3[3][3] = { { K_l3 * dyn2N, 0, 0 }, { 0, K_t3* dyn2N, 0 }, { 0, 0, K_t3* dyn2N } }; + + // c2 = 2/sqrt(6); + // c1 = 1/sqrt(6); + // c0 = 1/sqrt(2); + // c3 = 1/sqrt(3); + + double Phi10[3][3] = { { K_l3 * dyn2N, 0, 0 }, { 0, K_t3* dyn2N, 0 }, { 0, 0, K_t3* dyn2N } }; // Phi_nn3 + double Phi11[3][3]; + double Phi12[3][3]; + + // 4th neighbour + double r_j13[14][3] = { { -a_latt / sqrt3, 0, 0 }, + { a_latt * 0.5 / sqrt3, -a_latt * 0.5, 0 }, + { a_latt * 0.5 / sqrt3, a_latt * 0.5, 0 }, + { 0, a_latt, 0 }, + { -a_latt * sqrt3 / 2, a_latt / 2, 0 }, + { -a_latt * sqrt3 / 2, -a_latt / 2, 0 }, + { 0, -a_latt, 0 }, + { a_latt * sqrt3 / 2, -a_latt / 2, 0 }, + { a_latt * sqrt3 / 2, a_latt / 2, 0 }, + { 2 * a_latt / sqrt3, 0, 0 }, + { -a_latt / sqrt3, a_latt, 0 }, + { -a_latt / sqrt3, -a_latt, 0 }, + { 0, 0, c_latt / 2 }, + { 0, 0, -c_latt / 2 } }; // r_j12 + r_add3 Sublattice B and D do not have this coupling + + double Phi_nn4[3][3] = { { K_t4 * dyn2N, 0, 0 }, { 0, K_t4* dyn2N, 0 }, { 0, 0, K_l4* dyn2N } }; + + double Phi13[3][3] = { { K_t4 * dyn2N, 0, 0 }, { 0, K_t4* dyn2N, 0 }, { 0, 0, K_l4* dyn2N } }; // Phi_nn4; + double Phi14[3][3] = { { K_t4 * dyn2N, 0, 0 }, { 0, K_t4* dyn2N, 0 }, { 0, 0, K_l4* dyn2N } }; // Phi_nn4; + + // ----------------- END SPECIFIC FOR THE PG COPONENT -------------------------- + + /* TODO: Det er meget langsomt at bruge pointere. Erstat med flad memory. */ + // m(12*i+j) = M(i,j) + // 3*3 matrices multiplication + void + matrixmultiplication (double matrix1[3][3], double matrix2[3][3], double result[3][3]) { + // multiplies matrix1 by matrix2 and places the result in result + // double matrix1[3][3]={{*(array1), *(array1+1), *(array1+2)}, {*(array1+3), *(array1+4), *(array1+5)}, {*(array1+6), *(array1+7), *(array1+8)}}; + // double matrix2[3][3]={{*(array2), *(array2+1), *(array2+2)}, {*(array2+3), *(array2+4), *(array2+5)}, {*(array2+6), *(array2+7), *(array2+8)}}; + + result[0][0] = matrix2[0][0] * matrix1[0][0] + matrix2[1][0] * matrix1[0][1] + matrix2[2][0] * matrix1[0][2]; + result[0][1] = matrix2[0][1] * matrix1[0][0] + matrix2[1][1] * matrix1[0][1] + matrix2[2][1] * matrix1[0][2]; + result[0][2] = matrix2[0][2] * matrix1[0][0] + matrix2[1][2] * matrix1[0][1] + matrix2[2][2] * matrix1[0][2]; + result[1][0] = matrix2[0][0] * matrix1[1][0] + matrix2[1][0] * matrix1[1][1] + matrix2[2][0] * matrix1[1][2]; + result[1][1] = matrix2[0][1] * matrix1[1][0] + matrix2[1][1] * matrix1[1][1] + matrix2[2][1] * matrix1[1][2]; + result[1][2] = matrix2[0][2] * matrix1[1][0] + matrix2[1][2] * matrix1[1][1] + matrix2[2][2] * matrix1[1][2]; + result[2][0] = matrix2[0][0] * matrix1[2][0] + matrix2[1][0] * matrix1[2][1] + matrix2[2][0] * matrix1[2][2]; + result[2][1] = matrix2[0][1] * matrix1[2][0] + matrix2[1][1] * matrix1[2][1] + matrix2[2][1] * matrix1[2][2]; + result[2][2] = matrix2[0][2] * matrix1[2][0] + matrix2[1][2] * matrix1[2][1] + matrix2[2][2] * matrix1[2][2]; + return; + } + // TODO: (James Avery) Benchmark against pointer versions: matrix3mupltiplication or matrixmultiplication + // m(12*i+j) = M(i,j) + void + matrix3x3_multiply3 (const double A[9], const double B[9], const double C[9], double result[9]) { + // R_{i,j} = \sum_{k,l=1}^3 A_{i,k}*B_{k,l}*C_{l,j} + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) { + double ij_sum = 0; + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) + ij_sum += A[i * 3 + k] * B[k * 3 + l] * C[l * 3 + j]; + result[i * 3 + j] = ij_sum; + } + } -/* TODO: Det er meget langsomt at bruge pointere. Erstat med flad memory. */ -// m(12*i+j) = M(i,j) -//3*3 matrices multiplication -void matrixmultiplication( double matrix1[3][3], double matrix2[3][3], double result[3][3]) - { - // multiplies matrix1 by matrix2 and places the result in result -// double matrix1[3][3]={{*(array1), *(array1+1), *(array1+2)}, {*(array1+3), *(array1+4), *(array1+5)}, {*(array1+6), *(array1+7), *(array1+8)}}; -// double matrix2[3][3]={{*(array2), *(array2+1), *(array2+2)}, {*(array2+3), *(array2+4), *(array2+5)}, {*(array2+6), *(array2+7), *(array2+8)}}; - - result[0][0]=matrix2[0][0]*matrix1[0][0]+matrix2[1][0]*matrix1[0][1]+matrix2[2][0]*matrix1[0][2]; - result[0][1]=matrix2[0][1]*matrix1[0][0]+matrix2[1][1]*matrix1[0][1]+matrix2[2][1]*matrix1[0][2]; - result[0][2]=matrix2[0][2]*matrix1[0][0]+matrix2[1][2]*matrix1[0][1]+matrix2[2][2]*matrix1[0][2]; - result[1][0]=matrix2[0][0]*matrix1[1][0]+matrix2[1][0]*matrix1[1][1]+matrix2[2][0]*matrix1[1][2]; - result[1][1]=matrix2[0][1]*matrix1[1][0]+matrix2[1][1]*matrix1[1][1]+matrix2[2][1]*matrix1[1][2]; - result[1][2]=matrix2[0][2]*matrix1[1][0]+matrix2[1][2]*matrix1[1][1]+matrix2[2][2]*matrix1[1][2]; - result[2][0]=matrix2[0][0]*matrix1[2][0]+matrix2[1][0]*matrix1[2][1]+matrix2[2][0]*matrix1[2][2]; - result[2][1]=matrix2[0][1]*matrix1[2][0]+matrix2[1][1]*matrix1[2][1]+matrix2[2][1]*matrix1[2][2]; - result[2][2]=matrix2[0][2]*matrix1[2][0]+matrix2[1][2]*matrix1[2][1]+matrix2[2][2]*matrix1[2][2]; - - return; - } - -// TODO: (James Avery) Benchmark against pointer versions: matrix3mupltiplication or matrixmultiplication -// m(12*i+j) = M(i,j) -void matrix3x3_multiply3(const double A[9], const double B[9], const double C[9], double result[9]) -{ - // R_{i,j} = \sum_{k,l=1}^3 A_{i,k}*B_{k,l}*C_{l,j} - for(int i=0;i<3;i++) - for(int j=0;j<3;j++){ - double ij_sum = 0; - for(int k=0;k<3;k++) - for(int l=0;l<3;l++) - ij_sum += A[i*3+k]*B[k*3+l]*C[l*3+j]; - result[i*3+j] = ij_sum; - } -} - -/* TODO: (James Avery) It is slow to use points. Replace with flat memory. */ + /* TODO: (James Avery) It is slow to use points. Replace with flat memory. */ // m(12*i+j) = M(i,j) // Counter argument: (Kim Lefmann) pointers are now only used in the initialization. In the TRACE code, memory is flat - void matrix3multiplication(double matrix1[3][3], double matrix2[3][3], double matrix3[3][3], double result[3][3]) - { - // Performs the multiplication (matrix1*matrix2)*matrix3 and placed the result in result -//#define TEST_MULTIPLY - double temp[3][3]; - matrixmultiplication(matrix1,matrix2,temp); - matrixmultiplication(temp,matrix3,result); -#ifdef TEST_MULTIPLY - printf("Matrix3multiply called with \n"); - printf("matrix1 = ("); - for (int i=0; i<3; i++) - { - printf("( %g %g %g), ",matrix1[i][0], matrix1[i][1], matrix1[i][2]); - } - printf(")\n Matrix2 = ("); - for (int i=0; i<3; i++) - { - printf("( %g %g %g), ",matrix2[i][0], matrix2[i][1], matrix2[i][2]); - } - printf(")\n Matrix3 = ("); - for (int i=0; i<3; i++) - { - printf("( %g %g %g), ",matrix3[i][0], matrix3[i][1], matrix3[i][2]); - } - printf(")\n Temp = ("); - for (int i=0; i<3; i++) - { - printf("( %g %g %g), ",temp[i][0], temp[i][1], temp[i][2]); - } - printf(")\n Result = ("); - for (int i=0; i<3; i++) - { - printf("( %g %g %g), ",result[i][0], result[i][1], result[i][2]); - } - printf(")\n"); -#endif // TEST_MULTIPLY - - return; + void + matrix3multiplication (double matrix1[3][3], double matrix2[3][3], double matrix3[3][3], double result[3][3]) { + // Performs the multiplication (matrix1*matrix2)*matrix3 and placed the result in result + // #define TEST_MULTIPLY + double temp[3][3]; + matrixmultiplication (matrix1, matrix2, temp); + matrixmultiplication (temp, matrix3, result); + #ifdef TEST_MULTIPLY + printf ("Matrix3multiply called with \n"); + printf ("matrix1 = ("); + for (int i = 0; i < 3; i++) { + printf ("( %g %g %g), ", matrix1[i][0], matrix1[i][1], matrix1[i][2]); + } + printf (")\n Matrix2 = ("); + for (int i = 0; i < 3; i++) { + printf ("( %g %g %g), ", matrix2[i][0], matrix2[i][1], matrix2[i][2]); + } + printf (")\n Matrix3 = ("); + for (int i = 0; i < 3; i++) { + printf ("( %g %g %g), ", matrix3[i][0], matrix3[i][1], matrix3[i][2]); + } + printf (")\n Temp = ("); + for (int i = 0; i < 3; i++) { + printf ("( %g %g %g), ", temp[i][0], temp[i][1], temp[i][2]); } - + printf (")\n Result = ("); + for (int i = 0; i < 3; i++) { + printf ("( %g %g %g), ", result[i][0], result[i][1], result[i][2]); + } + printf (")\n"); + #endif // TEST_MULTIPLY -/* TODO: (James Avery) Replace by standard complex.h operation */ -complex complexexp_qdotr (double* q, double* rj) -{ -// double qrjtest - double qrj = *(q)**(rj)+*(q+1)**(rj+1)+*(q+2)**(rj+2); -// printf(" q= (%g %g %g), r = (%g %g %g) qrj = %g \n",q[0],q[1],q[2],rj[0],rj[1],rj[2],qrj); -// return (cexp(qrj)); // KL comment 030125 WHY does this not run? - return (cos(qrj)+I*sin(qrj)); -} + return; + } + /* TODO: (James Avery) Replace by standard complex.h operation */ + complex + complexexp_qdotr (double* q, double* rj) { + // double qrjtest + double qrj = *(q) * *(rj) + *(q + 1) * *(rj + 1) + *(q + 2) * *(rj + 2); + // printf(" q= (%g %g %g), r = (%g %g %g) qrj = %g \n",q[0],q[1],q[2],rj[0],rj[1],rj[2],qrj); + // return (cexp(qrj)); // KL comment 030125 WHY does this not run? + return (cos (qrj) + I * sin (qrj)); + } -void initialize_omega_q() -{ + void + initialize_omega_q () { // Make the matrix algebra that finalizes the force constant set-up - matrix3multiplication(Rot120,Phi_nn1,Rot120rev,Phi2); //Phi2=Rot120*Phi_nn1*Rot120^(-1) - matrix3multiplication(Rot120rev,Phi_nn1,Rot120,Phi3); //Phi3=Rot120^{-1}*Phi_nn1*Rot120 - matrix3multiplication(Rot60,Phi_nn2,Rot60rev,Phi5); //Phi5=Rot60*Phi_nn2*Rot60^{-1} - matrix3multiplication(Rot120,Phi_nn2,Rot120rev,Phi6); //Phi6=Rot120*Phi_nn2*Rot120^(-1); - matrix3multiplication(Rot180,Phi_nn2,Rot180,Phi7); //Phi7=Rot180*Phi_nn2*Rot180; - matrix3multiplication(Rot120rev,Phi_nn2,Rot120,Phi8); //Phi8=Rot120^{-1}*Phi_nn2*Rot120; - matrix3multiplication(Rot60rev,Phi_nn2,Rot60,Phi9); //Phi9=Rot60^{-1}*Phi_nn2*Rot60; - matrix3multiplication(Rot120,Phi_nn3,Rot120rev,Phi11);//Phi11=Rot120*Phi_nn3*Rot120^(-1); - matrix3multiplication(Rot120rev,Phi_nn3,Rot120,Phi12);//Phi12=Rot120^{-1}*Phi_nn3*Rot120; + matrix3multiplication (Rot120, Phi_nn1, Rot120rev, Phi2); // Phi2=Rot120*Phi_nn1*Rot120^(-1) + matrix3multiplication (Rot120rev, Phi_nn1, Rot120, Phi3); // Phi3=Rot120^{-1}*Phi_nn1*Rot120 + matrix3multiplication (Rot60, Phi_nn2, Rot60rev, Phi5); // Phi5=Rot60*Phi_nn2*Rot60^{-1} + matrix3multiplication (Rot120, Phi_nn2, Rot120rev, Phi6); // Phi6=Rot120*Phi_nn2*Rot120^(-1); + matrix3multiplication (Rot180, Phi_nn2, Rot180, Phi7); // Phi7=Rot180*Phi_nn2*Rot180; + matrix3multiplication (Rot120rev, Phi_nn2, Rot120, Phi8); // Phi8=Rot120^{-1}*Phi_nn2*Rot120; + matrix3multiplication (Rot60rev, Phi_nn2, Rot60, Phi9); // Phi9=Rot60^{-1}*Phi_nn2*Rot60; + matrix3multiplication (Rot120, Phi_nn3, Rot120rev, Phi11); // Phi11=Rot120*Phi_nn3*Rot120^(-1); + matrix3multiplication (Rot120rev, Phi_nn3, Rot120, Phi12); // Phi12=Rot120^{-1}*Phi_nn3*Rot120; return; -} + } -//---------------------------- START CENTRAL FUNCTION OMEGA-Q ---------------------------- + //---------------------------- START CENTRAL FUNCTION OMEGA-Q ---------------------------- + + double + omega_q (complex* parms, complex* eigenvectorgood) { - double omega_q(complex* parms, complex* eigenvectorgood) { - // _______________ DETERMINE Q ____________________________ - - double vi, vf, vv_x, vv_y, vv_z, vi_x, vi_y, vi_z; - double q, qx, qy, qz, Jq, hw_phonon, hw_neutron; - double h,k,l,hk_inplane, hk_outofplane; - int disp; - - for(int ii=0;ii<8;ii++) if(isnan(creal(parms[ii])) || isnan(cimag(parms[ii]))){ - fprintf(stderr,"omega_q: Array parms contains a nan at position %d.\n",ii); - fprintf(stderr,"parms = ["); for(int i=0;i<8;i++) fprintf(stderr,"%g + i%g%s",creal(parms[i]), cimag(parms[i]), i+1<8?", ":"]\n"); - abort(); - } - - vf=parms[0]; - vi=parms[1]; - vv_x=parms[2]; - vv_y=parms[3]; - vv_z=parms[4]; - vi_x=parms[5]; - vi_y=parms[6]; - vi_z=parms[7]; - h = parms[9]; - k = parms[10]; - l = parms[11]; - disp = (int)parms[12]; - - qx=V2K*(vi_x-vf*vv_x); - qy=V2K*(vi_y-vf*vv_y); - qz=V2K*(vi_z-vf*vv_z); - - if (disp==1) /* calculate dispersion directly */ - { - printf("Dispersion calculation: "); - qx = h*astar; - qy = k*astar; //correct only for qy=0 ! - qz = l*cstar; + + double vi, vf, vv_x, vv_y, vv_z, vi_x, vi_y, vi_z; + double q, qx, qy, qz, Jq, hw_phonon, hw_neutron; + double h, k, l, hk_inplane, hk_outofplane; + int disp; + + for (int ii = 0; ii < 8; ii++) + if (isnan (creal (parms[ii])) || isnan (cimag (parms[ii]))) { + fprintf (stderr, "omega_q: Array parms contains a nan at position %d.\n", ii); + fprintf (stderr, "parms = ["); + for (int i = 0; i < 8; i++) + fprintf (stderr, "%g + i%g%s", creal (parms[i]), cimag (parms[i]), i + 1 < 8 ? ", " : "]\n"); + abort (); } - - double qvec[3]={qx,qy,qz}; - - if (verbose>=4) - printf("Enter omega_q; q = (%g, %g, %g) \n",qx,qy,qz); - + + vf = parms[0]; + vi = parms[1]; + vv_x = parms[2]; + vv_y = parms[3]; + vv_z = parms[4]; + vi_x = parms[5]; + vi_y = parms[6]; + vi_z = parms[7]; + h = parms[9]; + k = parms[10]; + l = parms[11]; + disp = (int)parms[12]; + + qx = V2K * (vi_x - vf * vv_x); + qy = V2K * (vi_y - vf * vv_y); + qz = V2K * (vi_z - vf * vv_z); + + if (disp == 1) /* calculate dispersion directly */ + { + printf ("Dispersion calculation: "); + qx = h * astar; + qy = k * astar; // correct only for qy=0 ! + qz = l * cstar; + } + + double qvec[3] = { qx, qy, qz }; + + if (verbose >= 4) + printf ("Enter omega_q; q = (%g, %g, %g) \n", qx, qy, qz); + // ________________ END DETERMINE Q ________________- - - complex Phi_diag[3*3]; + complex Phi_diag[3 * 3]; int i; int j; for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - Phi_diag[i*3+j] = Phi1[i][j]+Phi2[i][j]+Phi3[i][j]+Phi4[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[3][0]))+Phi5[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[4][0]))+Phi6[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[5][0]))+Phi7[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[6][0]))+Phi8[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[7][0]))+Phi9[i][j]*(1-(*complexexp_qdotr)(&qvec[0], &r_j13[8][0]))+Phi10[i][j]+Phi11[i][j]+Phi12[i][j]; - } + for (j = 0; j < 3; j++) { + Phi_diag[i * 3 + j] = Phi1[i][j] + Phi2[i][j] + Phi3[i][j] + Phi4[i][j] * (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[3][0])) + + Phi5[i][j] * (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[4][0])) + Phi6[i][j] * (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[5][0])) + + Phi7[i][j] * (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[6][0])) + Phi8[i][j] * (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[7][0])) + + Phi9[i][j] * (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[8][0])) + Phi10[i][j] + Phi11[i][j] + Phi12[i][j]; + } } - - complex Phi_diag1[3*3]; + + complex Phi_diag1[3 * 3]; for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - Phi_diag1[i*3+j] = Phi13[i][j]+Phi14[i][j]; - } + for (j = 0; j < 3; j++) { + Phi_diag1[i * 3 + j] = Phi13[i][j] + Phi14[i][j]; + } } - - complex Phi_offdiag1[3*3]; - complex Phi_offdiag2[3*3]; + + complex Phi_offdiag1[3 * 3]; + complex Phi_offdiag2[3 * 3]; for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - Phi_offdiag1[i*3+j] =-Phi1[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[0][0]))-Phi2[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[1][0]))-Phi3[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[2][0]))-Phi10[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[9][0]))-Phi11[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[10][0]))-Phi12[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[11][0])); - Phi_offdiag2[i*3+j] = conj(Phi_offdiag1[i*3+j]); - } + for (j = 0; j < 3; j++) { + Phi_offdiag1[i * 3 + j] = -Phi1[i][j] * ((*complexexp_qdotr) (&qvec[0], &r_j13[0][0])) - Phi2[i][j] * ((*complexexp_qdotr) (&qvec[0], &r_j13[1][0])) + - Phi3[i][j] * ((*complexexp_qdotr) (&qvec[0], &r_j13[2][0])) - Phi10[i][j] * ((*complexexp_qdotr) (&qvec[0], &r_j13[9][0])) + - Phi11[i][j] * ((*complexexp_qdotr) (&qvec[0], &r_j13[10][0])) - Phi12[i][j] * ((*complexexp_qdotr) (&qvec[0], &r_j13[11][0])); + Phi_offdiag2[i * 3 + j] = conj (Phi_offdiag1[i * 3 + j]); + } } - - complex Phi_6Dim[6*6]; + + complex Phi_6Dim[6 * 6]; for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - Phi_6Dim[i*6+j] = Phi_diag[i*3+j]+Phi_diag1[i*3+j]; - Phi_6Dim[i*6+j+3] = Phi_offdiag1[i*3+j]; - Phi_6Dim[(i+3)*6+j] = Phi_offdiag2[i*3+j]; - Phi_6Dim[(i+3)*6+j+3] = Phi_diag[i*3+j]; // Because Phi_diag1 does not appear for the B,C, sublattices - } + for (j = 0; j < 3; j++) { + Phi_6Dim[i * 6 + j] = Phi_diag[i * 3 + j] + Phi_diag1[i * 3 + j]; + Phi_6Dim[i * 6 + j + 3] = Phi_offdiag1[i * 3 + j]; + Phi_6Dim[(i + 3) * 6 + j] = Phi_offdiag2[i * 3 + j]; + Phi_6Dim[(i + 3) * 6 + j + 3] = Phi_diag[i * 3 + j]; // Because Phi_diag1 does not appear for the B,C, sublattices + } } - - complex Phi_AC[3*3]; + + complex Phi_AC[3 * 3]; for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - Phi_AC[i*3+j] = -Phi13[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[12][0]))-Phi14[i][j]*((*complexexp_qdotr)(&qvec[0], &r_j13[13][0])); - } + for (j = 0; j < 3; j++) { + Phi_AC[i * 3 + j] = -Phi13[i][j] * ((*complexexp_qdotr) (&qvec[0], &r_j13[12][0])) - Phi14[i][j] * ((*complexexp_qdotr) (&qvec[0], &r_j13[13][0])); + } } - -// complex Zero_3Dim[3][3] = {{0 , 0 , 0},{0 , 0 , 0},{0 , 0 , 0}}; - complex Phi_6Dim_offdiag[6*6]; + // complex Zero_3Dim[3][3] = {{0 , 0 , 0},{0 , 0 , 0},{0 , 0 , 0}}; + + complex Phi_6Dim_offdiag[6 * 6]; for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - Phi_6Dim_offdiag[i*6+j] = Phi_AC[i*3+j]; - Phi_6Dim_offdiag[i*6+j+3] = (complex)0.0; - Phi_6Dim_offdiag[(i+3)*6+j] = (complex)0.0; - Phi_6Dim_offdiag[(i+3)*6+j+3] = (complex)0.0; - } + for (j = 0; j < 3; j++) { + Phi_6Dim_offdiag[i * 6 + j] = Phi_AC[i * 3 + j]; + Phi_6Dim_offdiag[i * 6 + j + 3] = (complex)0.0; + Phi_6Dim_offdiag[(i + 3) * 6 + j] = (complex)0.0; + Phi_6Dim_offdiag[(i + 3) * 6 + j + 3] = (complex)0.0; + } } - complex Phi_12Dim[12*12]; + complex Phi_12Dim[12 * 12]; for (i = 0; i < 6; i++) { - for (j = 0; j < 6; j++) { - Phi_12Dim[i*DIM+j] = Phi_6Dim[i*6+j]; - Phi_12Dim[i*DIM+j+6] = Phi_6Dim_offdiag[i*6+j]; - Phi_12Dim[(i+6)*DIM+j] = conj(Phi_6Dim_offdiag[i*6+j]); - Phi_12Dim[(i+6)*DIM+j+6] = Phi_6Dim[i*6+j]; - } + for (j = 0; j < 6; j++) { + Phi_12Dim[i * DIM + j] = Phi_6Dim[i * 6 + j]; + Phi_12Dim[i * DIM + j + 6] = Phi_6Dim_offdiag[i * 6 + j]; + Phi_12Dim[(i + 6) * DIM + j] = conj (Phi_6Dim_offdiag[i * 6 + j]); + Phi_12Dim[(i + 6) * DIM + j + 6] = Phi_6Dim[i * 6 + j]; + } } - -#ifdef TEST_MATRICES - printf("Phi_12Dim=\n"); + + #ifdef TEST_MATRICES + printf ("Phi_12Dim=\n"); for (i = 0; i < 12; i++) { - for (j = 0; j< 12; j++) { - if (j == 0) { - printf("{"); - } - printf("%g+(%g)i ", creal(Phi_12Dim[i*DIM+j]) , cimag((Phi_12Dim[i*DIM+j]))); - if (j == 11) { - printf("}\n"); - } + for (j = 0; j < 12; j++) { + if (j == 0) { + printf ("{"); + } + printf ("%g+(%g)i ", creal (Phi_12Dim[i * DIM + j]), cimag ((Phi_12Dim[i * DIM + j]))); + if (j == 11) { + printf ("}\n"); + } } } - - printf("Phi_6Dim=\n"); + + printf ("Phi_6Dim=\n"); for (i = 0; i < 6; i++) { - for (j = 0; j< 6; j++) { - if (j == 0) { - printf("{"); - } - printf("%g+(%g)i ", creal(Phi_6Dim[i*6+j]) , cimag((Phi_6Dim[i*6+j]))); - if (j == 5) { - printf("}\n"); - } + for (j = 0; j < 6; j++) { + if (j == 0) { + printf ("{"); + } + printf ("%g+(%g)i ", creal (Phi_6Dim[i * 6 + j]), cimag ((Phi_6Dim[i * 6 + j]))); + if (j == 5) { + printf ("}\n"); + } } } - printf("Phi_diag=\n"); + printf ("Phi_diag=\n"); for (i = 0; i < 3; i++) { - for (j = 0; j< 3; j++) { - if (j == 0) { - printf("{"); - } - printf("%g+(%g)i ", creal(Phi_diag[i*3+j]) , cimag((Phi_diag[i*3+j]))); - if (j == 2) { - printf("}\n"); - } + for (j = 0; j < 3; j++) { + if (j == 0) { + printf ("{"); + } + printf ("%g+(%g)i ", creal (Phi_diag[i * 3 + j]), cimag ((Phi_diag[i * 3 + j]))); + if (j == 2) { + printf ("}\n"); + } } } - printf("Phi2={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}}\n",Phi2[0][0],Phi2[0][1],Phi2[0][2],Phi2[1][0],Phi2[1][1],Phi2[1][2],Phi2[2][0],Phi2[2][1],Phi2[2][2]); - printf("Phi3={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi3[0][0],Phi3[0][1],Phi3[0][2],Phi3[1][0],Phi3[1][1],Phi3[1][2],Phi3[2][0],Phi3[2][1],Phi3[2][2]); - printf("Phi5={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi5[0][0],Phi5[0][1],Phi5[0][2],Phi5[1][0],Phi5[1][1],Phi5[1][2],Phi5[2][0],Phi5[2][1],Phi5[2][2]); - printf("Phi6={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi6[0][0],Phi6[0][1],Phi6[0][2],Phi6[1][0],Phi6[1][1],Phi6[1][2],Phi6[2][0],Phi6[2][1],Phi6[2][2]); - printf("Phi7={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi7[0][0],Phi7[0][1],Phi7[0][2],Phi7[1][0],Phi7[1][1],Phi7[1][2],Phi7[2][0],Phi7[2][1],Phi7[2][2]); - printf("Phi8={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi8[0][0],Phi8[0][1],Phi8[0][2],Phi8[1][0],Phi8[1][1],Phi8[1][2],Phi8[2][0],Phi8[2][1],Phi8[2][2]); - printf("Phi9={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi9[0][0],Phi9[0][1],Phi9[0][2],Phi9[1][0],Phi9[1][1],Phi9[1][2],Phi9[2][0],Phi9[2][1],Phi9[2][2]); - printf("Phi11={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi11[0][0],Phi11[0][1],Phi11[0][2],Phi11[1][0],Phi11[1][1],Phi11[1][2],Phi11[2][0],Phi11[2][1],Phi11[2][2]); - printf("Phi12={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi12[0][0],Phi12[0][1],Phi12[0][2],Phi12[1][0],Phi12[1][1],Phi12[1][2],Phi12[2][0],Phi12[2][1],Phi12[2][2]); - printf("Phi14={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}",Phi14[0][0],Phi14[0][1],Phi14[0][2],Phi14[1][0],Phi14[1][1],Phi14[1][2],Phi14[2][0],Phi14[2][1],Phi14[2][2]); - - printf(" 1-exp(q.r_j13[3]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[3][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[3][0]))); - printf(" 1-exp(q.r_j13[4]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[4][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[4][0]))); - printf(" 1-exp(q.r_j13[5]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[5][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[5][0]))); - printf(" 1-exp(q.r_j13[6]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[6][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[6][0]))); - printf(" 1-exp(q.r_j13[7]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[7][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[7][0]))); - printf(" 1-exp(q.r_j13[8]) = %lf + i %lf \n",creal(1-(*complexexp_qdotr)(&qvec[0], &r_j13[8][0])),cimag(1-(*complexexp_qdotr)(&qvec[0], &r_j13[8][0]))); - -#endif // TEST_MATRICES - - double eigenvalue[DIM]; - complex Q[DIM*DIM]; - complex matrix[DIM*DIM]; - int index[DIM]; - - - for (i=0; i= 6) - { - printf("Before bubblesort\n"); - for (i=0; i= 6) { + printf ("Before bubblesort\n"); + for (i = 0; i < DIM; i++) { + printf ("i = %i eigenvalue[i] = %g index[i] %i \n", i, eigenvalue[i], index[i]); } - } + } - // Sort eigenvalues from lowest to highest - bubblesort(DIM, eigenvalue, index); //mergesort(eigenvalue,DIM,index); - if (verbose >= 5) - { - printf("After bubblesort\n"); - for (i=0; i= 5) { + printf ("After bubblesort\n"); + for (i = 0; i < DIM; i++) { + printf ("i = %i eigenvalue[i] = %g index[i] %i \n", i, eigenvalue[i], index[i]); } - } + } + + double eigenvaluegood[DIM]; - double eigenvaluegood[DIM]; - - for (j = 0; j < DIM ; j++){ - eigenvectorgood[j] = Q[index[mode]*DIM + j]; /* Eigenvectors are columns of Q */ // No, they are rows! + for (j = 0; j < DIM; j++) { + eigenvectorgood[j] = Q[index[mode] * DIM + j]; /* Eigenvectors are columns of Q */ // No, they are rows! } - if (verbose>=7) - { - printf("Q = ["); - for (i = 0; i < DIM ; i++) - { - printf("\n ("); - for (j=0; j= 7) { + printf ("Q = ["); + for (i = 0; i < DIM; i++) { + printf ("\n ("); + for (j = 0; j < DIM; j++) + printf (" %g +i %g, ", creal (Q[i * DIM + j]), cimag (Q[i * DIM + j])); + printf (")"); } - - for (i=0; i=5) - { - printf("Diagonalization done, eigenvalues (energies) are: \n ("); - for (i=0; i= 4){ - printf("Return from omega_q \n"); - printf("hw_neutron = %g (vi = %g, vf = %g)\n",hw_neutron, vi, vf); - printf("hw_phonon = %g = eigenvaluegood[%d] = %g = eigenvalues[%d] = %g\n", - hw_phonon,mode,eigenvaluegood[mode],index[mode], eigenvalue[index[mode]]); -// printf("old parms = ["); for(int i=0;i<8;i++) fprintf(stderr,"%g+i%g%s",creal(parms[i]), cimag(parms[i]), i+1<8?", ":"]\n"); - } - - parms[8] = hw_phonon; - - if (disp==1) { - printf("mode = %d, (h,k,l) = (%g,%g,%g), hw_q = %g",mode,h,k,l,hw_phonon); - printf(" polarization: ( "); + + for (i = 0; i < DIM; i++) + eigenvaluegood[i] = sqrt (creal (eigenvalue[index[i]])) / sqrt (M) / (2 * PI * 1E12) * THz2meV; // convert to units of THz + + if (verbose >= 5) { + printf ("Diagonalization done, eigenvalues (energies) are: \n ("); + for (i = 0; i < DIM; i++) + printf (" %g, ", eigenvaluegood[i]); + printf (" )\n"); + printf ("mode = %i , eigenvectorgood[mode] = (", mode); for (j = 0; j < DIM; j++) - { - printf("%g + i %g ,",creal(eigenvectorgood[j]),cimag(eigenvectorgood[j])); + printf (" %g +i %g, ", creal (eigenvectorgood[j]), cimag (eigenvectorgood[j])); + printf (" )\n"); + } + + // ___________________ RETURN RESULTS TO CALLING FUNCTION __________________ + + hw_neutron = fabs (VS2E * (vi * vi - vf * vf)); // neutron energy transfer + // fabs used to find both positive and negative solutions, controlled by findroots() + hw_phonon = eigenvaluegood[mode]; // phonon energy + + if (verbose >= 4) { + printf ("Return from omega_q \n"); + printf ("hw_neutron = %g (vi = %g, vf = %g)\n", hw_neutron, vi, vf); + printf ("hw_phonon = %g = eigenvaluegood[%d] = %g = eigenvalues[%d] = %g\n", hw_phonon, mode, eigenvaluegood[mode], index[mode], eigenvalue[index[mode]]); + // printf("old parms = ["); for(int i=0;i<8;i++) fprintf(stderr,"%g+i%g%s",creal(parms[i]), cimag(parms[i]), i+1<8?", ":"]\n"); + } + + parms[8] = hw_phonon; + + if (disp == 1) { + printf ("mode = %d, (h,k,l) = (%g,%g,%g), hw_q = %g", mode, h, k, l, hw_phonon); + printf (" polarization: ( "); + for (j = 0; j < DIM; j++) { + printf ("%g + i %g ,", creal (eigenvectorgood[j]), cimag (eigenvectorgood[j])); } - printf(")\n"); - } - if (disp==2) { - hk_inplane = qx/astar; - hk_outofplane = qy/astar; // Derivation assume astar along x and 60 degrees between astar and bstar: + printf (")\n"); + } + if (disp == 2) { + hk_inplane = qx / astar; + hk_outofplane = qy / astar; // Derivation assume astar along x and 60 degrees between astar and bstar: // h_ip = h + k/2 h_op = sqrt(3)*k/2 => k = 2/sqrt(3)*h_op h = h_ip - h_op/sqrt(3) - k = 2/sqrt(3)*hk_outofplane; - h = hk_inplane - hk_outofplane/sqrt(3); - l = qz/cstar; - printf("mode = %d, vf = %g, (h,k,l) = (%g,%g,%g), hw_q = %g",mode,vf,h,k,l,hw_phonon); - printf(" polarization: (%g + i %g , %g + i %g , %g + i %g) \n",creal(eigenvectorgood[0]),cimag(eigenvectorgood[0]),creal(eigenvectorgood[1]),cimag(eigenvectorgood[1]), creal(eigenvectorgood[2]),cimag(eigenvectorgood[2])); - } + k = 2 / sqrt (3) * hk_outofplane; + h = hk_inplane - hk_outofplane / sqrt (3); + l = qz / cstar; + printf ("mode = %d, vf = %g, (h,k,l) = (%g,%g,%g), hw_q = %g", mode, vf, h, k, l, hw_phonon); + printf (" polarization: (%g + i %g , %g + i %g , %g + i %g) \n", creal (eigenvectorgood[0]), cimag (eigenvectorgood[0]), creal (eigenvectorgood[1]), + cimag (eigenvectorgood[1]), creal (eigenvectorgood[2]), cimag (eigenvectorgood[2])); + } - return (hw_phonon - hw_neutron); -} + return (hw_phonon - hw_neutron); + } -//--------------------- END CENTRAL FUNCTION OMEGA-Q ------------------------ + //--------------------- END CENTRAL FUNCTION OMEGA-Q ------------------------ -// ------------------START OF THE GENERAL ROOT-FINDING CODE-------------------- + // ------------------START OF THE GENERAL ROOT-FINDING CODE-------------------- -double last_fh, last_x2; -double zridd(double (*func)(complex*,complex*), double x1, double x2, complex *parms, complex* vector, double xacc) - { - int j; - double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; + double last_fh, last_x2; + double + zridd (double (*func) (complex*, complex*), double x1, double x2, complex* parms, complex* vector, double xacc) { + int j; + double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; - // fprintf(stderr,"i zridd(%g,%g, parms,%g)\n",x1,x2,xacc); - parms[0]=x1; - // fprintf(stderr,"zridd fl = omega_q\n"); - if (xacc>0) - { - fl=(*func)(parms,vector); - } - else - { - fl = last_fh; - xacc = -xacc; - if (fabs(x1-last_x2)>xacc) // KL comment 030125 WHY does the library routine zridd use complex abs() function ?? changed to fabs() - printf("*** error in zridd *** x1: %g last_x2: %g \n",x1,last_x2); + // fprintf(stderr,"i zridd(%g,%g, parms,%g)\n",x1,x2,xacc); + parms[0] = x1; + // fprintf(stderr,"zridd fl = omega_q\n"); + if (xacc > 0) { + fl = (*func) (parms, vector); + } else { + fl = last_fh; + xacc = -xacc; + if (fabs (x1 - last_x2) > xacc) // KL comment 030125 WHY does the library routine zridd use complex abs() function ?? changed to fabs() + printf ("*** error in zridd *** x1: %g last_x2: %g \n", x1, last_x2); + } + // fprintf(stderr,"fl = %g\n",fl); + parms[0] = x2; + // fprintf(stderr,"zridd fh = omega_q; parms[0] = x2 = %g+i%g = %g\n",creal(parms[0]), cimag(parms[0]),x2); + last_fh = fh = (*func) (parms, vector); + last_x2 = x2; + // fprintf(stderr,"fh = %g\n",fh); + if (fl * fh >= 0) { + if (fl == 0) + return x1; + if (fh == 0) + return x2; + return UNUSED; + } else { + xl = x1; + xh = x2; + // printf("zridd sign change: v_low : %g v_high: %g \n",xl,xh); + ans = UNUSED; + for (j = 1; j < MAXRIDD; j++) { + xm = 0.5 * (xl + xh); + parms[0] = xm; + // fprintf(stderr,"zridd fm = omega_q\n"); + fm = (*func) (parms, vector); + s = sqrt (fm * fm - fl * fh); + if (s == 0.0) + return ans; + xnew = xm + (xm - xl) * ((fl >= fh ? 1.0 : -1.0) * fm / s); + if (fabs (xnew - ans) <= xacc) + return ans; + ans = xnew; + parms[0] = ans; + // fprintf(stderr,"s = %g; fl, fm, fh = %g,%g,%g; fnew = %g; fm*fm-fl*fh = %g\n", + // s, fl, fm,fh, fnew, fm*fm-fl*fh); + // fprintf(stderr,"zridd fnew = omega_q\n"); + fnew = (*func) (parms, vector); + if (fnew == 0.0) + return ans; + if (fabs (fm) * SIGN (fnew) != fm) { + xl = xm; + fl = fm; + xh = ans; + fh = fnew; + } else if (fabs (fl) * SIGN (fnew) != fl) { + xh = ans; + fh = fnew; + } else if (fabs (fh) * SIGN (fnew) != fh) { + xl = ans; + fl = fnew; + } else + fatalerror ("never get here in zridd"); + if (fabs (xh - xl) <= xacc) + return ans; } - // fprintf(stderr,"fl = %g\n",fl); - parms[0]=x2; - // fprintf(stderr,"zridd fh = omega_q; parms[0] = x2 = %g+i%g = %g\n",creal(parms[0]), cimag(parms[0]),x2); - last_fh=fh=(*func)(parms,vector); - last_x2 = x2; - // fprintf(stderr,"fh = %g\n",fh); - if (fl*fh >= 0) - { - if (fl==0) return x1; - if (fh==0) return x2; - return UNUSED; + fatalerror ("zridd exceeded maximum iterations"); + } + return 0.0; /* Never get here */ + } + + #pragma acc routine + double + zridd_gpu (double x1, double x2, complex* parms, complex* vector, double xacc) { + int j; + double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; + + parms[0] = x1; + // fprintf(stderr,"fl=omega_q(%g+i%g)\n",creal(parms[0]),cimag(parms[0])); + fl = omega_q (parms, vector); + parms[0] = x2; + // fprintf(stderr,"fh=omega_q(%g+i%g)\n",creal(parms[0]),cimag(parms[0])); + fh = omega_q (parms, vector); + if (fl * fh >= 0) { + if (fl == 0) + return x1; + if (fh == 0) + return x2; + return UNUSED; + } else { + xl = x1; + xh = x2; + ans = UNUSED; + for (j = 1; j < MAXRIDD; j++) { + xm = 0.5 * (xl + xh); + parms[0] = xm; + // fprintf(stderr,"fm=omega_q(%g+i%g)\n",creal(parms[0]),cimag(parms[0])); + fm = omega_q (parms, vector); + s = sqrt (fm * fm - fl * fh); + if (s == 0.0) + return ans; + xnew = xm + (xm - xl) * ((fl >= fh ? 1.0 : -1.0) * fm / s); + if (fabs (xnew - ans) <= xacc) + return ans; + ans = xnew; + parms[0] = ans; + // fprintf(stderr,"fnew=omega_q(%g+i%g)\n",creal(parms[0]),cimag(parms[0])); + fnew = omega_q (parms, vector); + if (fnew == 0.0) + return ans; + if (fabs (fm) * SIGN (fnew) != fm) { + xl = xm; + fl = fm; + xh = ans; + fh = fnew; + } else if (fabs (fl) * SIGN (fnew) != fl) { + xh = ans; + fh = fnew; + } else if (fabs (fh) * SIGN (fnew) != fh) { + xl = ans; + fl = fnew; + } else + fatalerror ("never get here in zridd"); + if (fabs (xh - xl) <= xacc) + return ans; } + fatalerror ("zridd exceeded maximum iterations"); + } + return 0.0; /* Never get here */ + } + + #define ROOTACC 1e-8 + + int + findroots (double brack_low, double brack_mid, double brack_high, double* list, int* index, double (*f) (complex*, complex*), complex* parms, complex* vector, + int low_steps, int high_steps) + // *f is here omega_q + // *parms is here p_call + { + double root, range; + int i; + + if (verbose == 2) + printf ("Enter findroots() \n"); + + // First, find roots in energy loss side, if any + range = brack_mid - brack_low; + for (i = 0; i < (low_steps - 1); i++) { + if (i == 0) + root = zridd (f, brack_low + range * i / low_steps, brack_low + range * (i + 1) / low_steps, (complex*)parms, (complex*)vector, ROOTACC); else - { - xl=x1; - xh=x2; -// printf("zridd sign change: v_low : %g v_high: %g \n",xl,xh); - ans=UNUSED; - for (j=1; j= fh ? 1.0 : -1.0)*fm/s); - if (fabs(xnew-ans) <= xacc) - return ans; - ans=xnew; - parms[0]=ans; - // fprintf(stderr,"s = %g; fl, fm, fh = %g,%g,%g; fnew = %g; fm*fm-fl*fh = %g\n", - // s, fl, fm,fh, fnew, fm*fm-fl*fh); - //fprintf(stderr,"zridd fnew = omega_q\n"); - fnew=(*func)(parms,vector); - if (fnew == 0.0) return ans; - if (fabs(fm)*SIGN(fnew) != fm) - { - xl=xm; - fl=fm; - xh=ans; - fh=fnew; - } - else - if (fabs(fl)*SIGN(fnew) != fl) - { - xh=ans; - fh=fnew; - } - else - if(fabs(fh)*SIGN(fnew) != fh) - { - xl=ans; - fl=fnew; - } - else - fatalerror("never get here in zridd"); - if (fabs(xh-xl) <= xacc) - return ans; - } - fatalerror("zridd exceeded maximum iterations"); + root = zridd (f, brack_low + range * i / low_steps, brack_low + range * (i + 1) / low_steps, (complex*)parms, (complex*)vector, + -ROOTACC); // Re-use the previous fh value + if (root != UNUSED) { + list[(*index)++] = root; + if (verbose >= 4) + printf ("--- findroots returned a root on the energy loss side: vf = %g \n", root); } - return 0.0; /* Never get here */ } -#pragma acc routine -double zridd_gpu(double x1, double x2, complex *parms, complex* vector, double xacc) + range = (brack_mid - brack_low) / (double)low_steps; + for (i = 0; i < low_steps; i++) // Small steps close to zero energy transfer { - int j; - double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; - - parms[0]=x1; - // fprintf(stderr,"fl=omega_q(%g+i%g)\n",creal(parms[0]),cimag(parms[0])); - fl=omega_q(parms,vector); - parms[0]=x2; - // fprintf(stderr,"fh=omega_q(%g+i%g)\n",creal(parms[0]),cimag(parms[0])); - fh=omega_q(parms,vector); - if (fl*fh >= 0) - { - if (fl==0) return x1; - if (fh==0) return x2; - return UNUSED; - } + if (low_steps == 1) // then the previous loop did not execute + root = zridd (f, brack_low + range * (low_steps - 1 + i / (double)low_steps), brack_low + range * (low_steps - 1 + (i + 1) / (double)low_steps), + (complex*)parms, (complex*)vector, ROOTACC); else - { - xl=x1; - xh=x2; - ans=UNUSED; - for (j=1; j= fh ? 1.0 : -1.0)*fm/s); - if (fabs(xnew-ans) <= xacc) - return ans; - ans=xnew; - parms[0]=ans; - // fprintf(stderr,"fnew=omega_q(%g+i%g)\n",creal(parms[0]),cimag(parms[0])); - fnew=omega_q(parms,vector); - if (fnew == 0.0) return ans; - if (fabs(fm)*SIGN(fnew) != fm) - { - xl=xm; - fl=fm; - xh=ans; - fh=fnew; - } - else - if (fabs(fl)*SIGN(fnew) != fl) - { - xh=ans; - fh=fnew; - } - else - if(fabs(fh)*SIGN(fnew) != fh) - { - xl=ans; - fl=fnew; - } - else - fatalerror("never get here in zridd"); - if (fabs(xh-xl) <= xacc) - return ans; - } - fatalerror("zridd exceeded maximum iterations"); + root = zridd (f, brack_low + range * (low_steps - 1 + i / (double)low_steps), brack_low + range * (low_steps - 1 + (i + 1) / (double)low_steps), + (complex*)parms, (complex*)vector, -ROOTACC); // Re-use the previous fh value + if (root != UNUSED) { + list[(*index)++] = root; + if (verbose >= 4) + printf ("--- findroots returned a root on the energy loss side: vf = %g \n", root); } - return 0.0; /* Never get here */ } - -#define ROOTACC 1e-8 - int findroots(double brack_low, double brack_mid, double brack_high, double *list, int* index, double (*f)(complex*,complex*), complex *parms, complex* vector, int low_steps, int high_steps) - // *f is here omega_q - // *parms is here p_call + // Second, find roots in energy gain side, there is always some + range = (brack_high - brack_mid) / (double)high_steps; + for (i = 0; i < high_steps; i++) // Close to zero: small steps { - double root,range; - int i; - - if (verbose==2) - printf("Enter findroots() \n"); - - // First, find roots in energy loss side, if any - range = brack_mid-brack_low; - for (i=0; i<(low_steps-1); i++) - { - if (i==0) - root = zridd(f, brack_low+range*i/low_steps, - brack_low+range*(i+1)/low_steps, - (complex *)parms, (complex *)vector, ROOTACC); - else - root = zridd(f, brack_low+range*i/low_steps, - brack_low+range*(i+1)/low_steps, - (complex *)parms, (complex *)vector, -ROOTACC); // Re-use the previous fh value - if (root != UNUSED) - { - list[(*index)++]=root; - if (verbose >=4) - printf("--- findroots returned a root on the energy loss side: vf = %g \n",root); - } - } - - range = (brack_mid-brack_low)/(double)low_steps; - for (i=0; i=4) - printf("--- findroots returned a root on the energy loss side: vf = %g \n",root); - } - } - - // Second, find roots in energy gain side, there is always some - range = (brack_high-brack_mid)/(double)high_steps; - for (i=0; i= 4) - printf("*** findroots returned a root on the energy gain side: vf = %g \n",root); + printf ("*** findroots returned a root on the energy gain side: vf = %g \n", root); } - } - - range = brack_high-brack_mid; - for (i=1; i= 4) - printf("*** findroots returned a root on the energy gain side: vf = %g \n",root); + printf ("*** findroots returned a root on the energy gain side: vf = %g \n", root); } - } - } - -#pragma acc routine - int findroots_gpu(double brack_low, double brack_mid, double brack_high, double *list, int* index, double *parms, complex* vector, int low_steps, int high_steps) - { - double root,range=brack_mid-brack_low; - int i; - - for (i=0; i 0 && yheight) shape=0; /* cylinder */ + int i, j; + if (xwidth && yheight && zdepth) + shape = 1; /* box */ + else if (radius > 0 && yheight) + shape = 0; /* cylinder */ - V_rho = 1/(a_latt*a_latt*c_latt*sqrt3/2); // (unit: Å^-3) + V_rho = 1 / (a_latt * a_latt * c_latt * sqrt3 / 2); // (unit: Å^-3) V_my_s = (V_rho * 100 * sigma_inc); V_my_a_v = (V_rho * 100 * sigma_abs * 2200); /* now compute target coords if a component index is supplied */ - if (!target_index && !target_x && !target_y && !target_z) target_index=1; - if (target_index){ + if (!target_index && !target_x && !target_y && !target_z) + target_index = 1; + if (target_index) { Coords ToTarget; - ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+target_index), POS_A_CURRENT_COMP); - ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget); - coords_get(ToTarget, &target_x, &target_y, &target_z); + ToTarget = coords_sub (POS_A_COMP_INDEX (INDEX_CURRENT_COMP + target_index), POS_A_CURRENT_COMP); + ToTarget = rot_apply (ROT_A_CURRENT_COMP, ToTarget); + coords_get (ToTarget, &target_x, &target_y, &target_z); } if (!(target_x || target_y || target_z)) { - printf("Phonon_BkV_PG: %s: The target is not defined. Using direct beam (Z-axis).\n", - NAME_CURRENT_COMP); - target_z=1; + printf ("Phonon_BkV_PG: %s: The target is not defined. Using direct beam (Z-axis).\n", NAME_CURRENT_COMP); + target_z = 1; } - initialize_omega_q(); + initialize_omega_q (); verbose = verbose_input; - // m(12*i+j) = M(i,j) - // OK for: Rot120 - for (i=0; i<3; i++) - for (j=0; j<3; j++) - { - printf("i,j: %i, %i ROT120: %g ROT120FLAT; %g \n",i,j,Rot120[i][j],Rot120_flat[3*i+j]); + // m(12*i+j) = M(i,j) + // OK for: Rot120 + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) { + printf ("i,j: %i, %i ROT120: %g ROT120FLAT; %g \n", i, j, Rot120[i][j], Rot120_flat[3 * i + j]); } %} TRACE %{ - double t0, t1, t2, t3; /* Entry/exit times for shape */ - double v_i, v_f; /* Neutron velocities: initial, final */ - double vx_i, vy_i, vz_i; /* Neutron initial velocity vector */ - double dt0, dt; /* Flight times through sample */ - double l_full; /* Flight path length for non-scattered neutron */ - double l_i, l_o; /* Flight path lenght in/out for scattered neutron */ - double my_a_i; /* Initial attenuation factor */ - double my_a_f; /* Final attenuation factor */ - double solid_angle; /* Solid angle of target as seen from scattering point */ - double aim_x=0, aim_y=0, aim_z=1; /* Position of target relative to scattering point */ - double omega; /* Energy transfer */ - double qsquare; /* Square of the scattering vector */ - double bose_factor; /* Calculated value of the Bose factor */ - int nf, index; /* Number of allowed final velocities */ - double vf_list[20]; /* List of allowed final velocities */ - double J_factor, J0; /* Jacobian from delta fnc.s in cross section; elastic J-factor */ - double f1, f2; /* Probed values of omega_q minus omega */ - double p_att,p_MC,p_ph1,p_ph2,p_ph3,p_J, p_tot; /* Temporary multipliers */ - complex Gn; /* Phonon structure factor */ - complex q_dot_e; /* q dot e */ - double q_dot_r; /* q dot r */ - double qh,qk,ql,qhkx,qhky,qh_Bragg,qk_Bragg,ql_Bragg,qh_res,qk_res,qh_add,qk_add; /* components of q */ - double dist_00, dist_01, dist_10, dist_11, dist_min; /* used to calculate nearest Bragg point */ - int i,j, intersect; - - if (shape == 0) - intersect = cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius, yheight); - else if (shape == 1) - intersect = box_intersect(&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth); - - if(intersect) - { - // fprintf(stderr,"intersect\n"); - if(t0 < 0) + double t0, t1, t2, t3; /* Entry/exit times for shape */ + double v_i, v_f; /* Neutron velocities: initial, final */ + double vx_i, vy_i, vz_i; /* Neutron initial velocity vector */ + double dt0, dt; /* Flight times through sample */ + double l_full; /* Flight path length for non-scattered neutron */ + double l_i, l_o; /* Flight path lenght in/out for scattered neutron */ + double my_a_i; /* Initial attenuation factor */ + double my_a_f; /* Final attenuation factor */ + double solid_angle; /* Solid angle of target as seen from scattering point */ + double aim_x = 0, aim_y = 0, aim_z = 1; /* Position of target relative to scattering point */ + double omega; /* Energy transfer */ + double qsquare; /* Square of the scattering vector */ + double bose_factor; /* Calculated value of the Bose factor */ + int nf, index; /* Number of allowed final velocities */ + double vf_list[20]; /* List of allowed final velocities */ + double J_factor, J0; /* Jacobian from delta fnc.s in cross section; elastic J-factor */ + double f1, f2; /* Probed values of omega_q minus omega */ + double p_att, p_MC, p_ph1, p_ph2, p_ph3, p_J, p_tot; /* Temporary multipliers */ + complex Gn; /* Phonon structure factor */ + complex q_dot_e; /* q dot e */ + double q_dot_r; /* q dot r */ + double qh, qk, ql, qhkx, qhky, qh_Bragg, qk_Bragg, ql_Bragg, qh_res, qk_res, qh_add, qk_add; /* components of q */ + double dist_00, dist_01, dist_10, dist_11, dist_min; /* used to calculate nearest Bragg point */ + int i, j, intersect; + + if (shape == 0) + intersect = cylinder_intersect (&t0, &t3, x, y, z, vx, vy, vz, radius, yheight); + else if (shape == 1) + intersect = box_intersect (&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth); + + if (intersect) { + // fprintf(stderr,"intersect\n"); + if (t0 < 0) ABSORB; /* Neutron came from the sample or begins inside */ - if (verbose>=2) - printf("neutron entered Phonon_BvK_PG \n"); + if (verbose >= 2) + printf ("neutron entered Phonon_BvK_PG \n"); /* Neutron enters at t=t0. */ - dt0 = t3-t0; /* Time in sample */ - v_i = sqrt(vx*vx + vy*vy + vz*vz); + dt0 = t3 - t0; /* Time in sample */ + v_i = sqrt (vx * vx + vy * vy + vz * vz); l_full = v_i * dt0; /* Length of path through sample if not scattered */ - dt = rand01()*dt0; /* Time of scattering (relative to t0) */ - l_i = v_i*dt; /* Penetration in sample at scattering */ - vx_i=vx; - vy_i=vy; - vz_i=vz; - PROP_DT(dt+t0); /* Point of scattering */ - - aim_x = target_x-x; /* Vector pointing at target (e.g. analyzer) */ - aim_y = target_y-y; - aim_z = target_z-z; - - if(focus_aw && focus_ah) { - randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle, - aim_x, aim_y, aim_z, focus_aw, focus_ah, ROT_A_CURRENT_COMP); - } else if(focus_xw && focus_yh) { - randvec_target_rect(&vx, &vy, &vz, &solid_angle, - aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP); + dt = rand01 () * dt0; /* Time of scattering (relative to t0) */ + l_i = v_i * dt; /* Penetration in sample at scattering */ + vx_i = vx; + vy_i = vy; + vz_i = vz; + PROP_DT (dt + t0); /* Point of scattering */ + + aim_x = target_x - x; /* Vector pointing at target (e.g. analyzer) */ + aim_y = target_y - y; + aim_z = target_z - z; + + if (focus_aw && focus_ah) { + randvec_target_rect_angular (&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_aw, focus_ah, ROT_A_CURRENT_COMP); + } else if (focus_xw && focus_yh) { + randvec_target_rect (&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP); } else { - randvec_target_sphere(&vx,&vy,&vz,&solid_angle,aim_x,aim_y,aim_z, focus_r); + randvec_target_sphere (&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_r); } - NORM(vx, vy, vz); - nf=0; + NORM (vx, vy, vz); + nf = 0; if (mode_input == 12) - mode = round(12*rand01()-0.5); // Select the phonon mode + mode = round (12 * rand01 () - 0.5); // Select the phonon mode else - mode = mode_input; // Use the given mode - if ((mode < 0) || (mode > 11)) // Undefined mode specification + mode = mode_input; // Use the given mode + if ((mode < 0) || (mode > 11)) // Undefined mode specification { - printf("mode = %d ",mode); - nrerror("illegal value of mode"); + printf ("mode = %d ", mode); + nrerror ("illegal value of mode"); } - p_call[0]=-1; // in the end this will be v_f - p_call[1]=v_i; - p_call[2]=vx; - p_call[3]=vy; - p_call[4]=vz; - p_call[5]=vx_i; - p_call[6]=vy_i; - p_call[7]=vz_i; - p_call[9]=hh; - p_call[10]=kk; - p_call[11]=ll; - p_call[12] = (double)dispersion; - - if (dispersion==1) - { - f1=omega_q(p_call,eigenvectormode); - p_call[12]=0.0; - } + p_call[0] = -1; // in the end this will be v_f + p_call[1] = v_i; + p_call[2] = vx; + p_call[3] = vy; + p_call[4] = vz; + p_call[5] = vx_i; + p_call[6] = vy_i; + p_call[7] = vz_i; + p_call[9] = hh; + p_call[10] = kk; + p_call[11] = ll; + p_call[12] = (double)dispersion; + + if (dispersion == 1) { + f1 = omega_q (p_call, eigenvectormode); + p_call[12] = 0.0; + } -/* if (dispersion==2) - { - for(i=0; i=2) - printf("Call findroots \n"); - findroots(0, v_i, v_i+V_HIGH, vf_list, &nf, omega_q, p_call, eigenvectormode, e_steps_low, e_steps_high); - if (verbose >=2) - printf( "Findroots returned %d roots \n",nf); + if (verbose >= 2) + printf ("Call findroots \n"); + findroots (0, v_i, v_i + V_HIGH, vf_list, &nf, omega_q, p_call, eigenvectormode, e_steps_low, e_steps_high); + if (verbose >= 2) + printf ("Findroots returned %d roots \n", nf); #else - findroots_gpu(0, v_i, v_i+V_HIGH, vf_list, &nf, p_call, eigenvectormode, e_steps_low, e_steps_high); + findroots_gpu (0, v_i, v_i + V_HIGH, vf_list, &nf, p_call, eigenvectormode, e_steps_low, e_steps_high); #endif - if (verbose>=3) - { - printf("Findroots done, mode %i , last phonon energy %g \n",mode,creal(p_call[8])); + if (verbose >= 3) { + printf ("Findroots done, mode %i , last phonon energy %g \n", mode, creal (p_call[8])); } - + // ________________ ALL FROM HERE IS CALCULATION OF INTENSITY - SHOULD BE REVISITED ___________________ - - index=(int)floor(rand01()*nf); // Select random solution - v_f=vf_list[index]; - p_call[0]=v_f; // transfer choice of v_f to omega_q - - f1=omega_q(p_call,eigenvectormode); - if (verbose >= 2) - { - printf("Chosen solution is %i; v_f = %g, hw = %g, energy match is %g; eigenvector is: (", index, v_f, creal(p_call[8]), f1); - for (j=0; j= 3) - { - printf("q transforms in r.l.u to (h, k, l) = (%g, %g, %g). Nearest Bragg point: (%g, %g, %g)\n",qh,qk,ql,qh_Bragg,qk_Bragg,ql_Bragg); - printf("--------------------------------------------\n"); - } + f1 = omega_q (p_call, eigenvectormode); + if (verbose >= 2) { + printf ("Chosen solution is %i; v_f = %g, hw = %g, energy match is %g; eigenvector is: (", index, v_f, creal (p_call[8]), f1); + for (j = 0; j < DIM; j++) + printf (" %g +i %g, ", creal (eigenvectormode[j]), cimag (eigenvectormode[j])); + printf (" )\n"); + printf ("--------------------------------------------\n"); + } + p_call[0] = v_f - DV; + // fprintf(stderr,"f1=omega_q(%g+i%g)\n",creal(p_call[0]),cimag(p_call[0])); + f1 = omega_q (p_call, eigenvectormode); + // fprintf(stderr,"retur fra f1=omega_q\n"); + p_call[0] = v_f + DV; + // fprintf(stderr,"f2=omega_q(%g+i%g)\n",creal(p_call[0]),cimag(p_call[0])); + f2 = omega_q (p_call, eigenvectormode); + // fprintf(stderr,"retur fra f2=omega_q\n"); + J_factor = 1.6022E-22 * 1E-10 * fabs (f2 - f1) / (2 * DV * V2K); // unit: meV Å converted to SI units + omega = VS2E * (v_i * v_i - v_f * v_f); // unit: meV + vx *= v_f; + vy *= v_f; + vz *= v_f; + + q_x = q[0] = V2K * (vx_i - vx); + q_y = q[1] = V2K * (vy_i - vy); + q_z = q[2] = V2K * (vz_i - vz); + + ql = q_z / cstar; + qhkx = q_x / astar; // These are the Cartesian coordinates + qhky = q_y / astar; // These are the Cartesian coordinates + qk = 2 * qhky / sqrt3; // transform to hexagonal coordinates + qh = qhkx - qhky / sqrt3; // transform to hexagonal coordinates + + // Now find out which Bragg peak is closest in hexagonal plane + + ql_Bragg = round (ql); // We know what to do along the c-axis + + qh_Bragg = floor (qh); // Divide into integer and residual part + qh_res = qh - qh_Bragg; + qk_Bragg = floor (qk); + qk_res = qk - qk_Bragg; + + dist_00 = qh_res * qh_res + qk_res * qk_res + qh_res * qk_res; // Square distance in hexagonal coordinates to (0,0) + dist_min = dist_00; + qh_add = 0; + qk_add = 0; + + dist_01 = qh_res * qh_res + (qk_res - 1) * (qk_res - 1) + qh_res * (qk_res - 1); // Square distance in hexagonal coordinates to (0,1) + if (dist_01 < dist_min) { + dist_min = dist_01; + qh_add = 0; + qk_add = 1; + } + dist_10 = (qh_res - 1) * (qh_res - 1) + qk_res * qk_res + (qh_res - 1) * qk_res; // Square distance in hexagonal coordinates to (1,0) + if (dist_10 < dist_min) { + dist_min = dist_10; + qh_add = 1; + qk_add = 0; + } + dist_11 = (qh_res - 1) * (qh_res - 1) + (qk_res - 1) * (qk_res - 1) + (qh_res - 1) * (qk_res - 1); // Square distance in hexagonal coordinates to (1,1) + if (dist_11 < dist_min) { + dist_min = dist_11; + qh_add = 1; + qk_add = 1; + } + qh_Bragg += qh_add; // Potentially correct to the right BZ + qk_Bragg += qk_add; - if (dispersion==1) - printf("Dispersion found: (h,k,l) = (%g, %g, %g), hw = %g \n",q[0]/astar,q[1]/astar,q[2]/cstar,omega); - - qsquare=1; - Gn = 0; + if (verbose >= 3) { + printf ("q transforms in r.l.u to (h, k, l) = (%g, %g, %g). Nearest Bragg point: (%g, %g, %g)\n", qh, qk, ql, qh_Bragg, qk_Bragg, ql_Bragg); + printf ("--------------------------------------------\n"); + } + + if (dispersion == 1) + printf ("Dispersion found: (h,k,l) = (%g, %g, %g), hw = %g \n", q[0] / astar, q[1] / astar, q[2] / cstar, omega); - double q_Bragg[3]; // Bragg point in Å-1; Cartesian coordinates - q_Bragg[0]=qh_Bragg*astar+qk_Bragg*astar/2.0; - q_Bragg[1]=qk_Bragg*astar*2.0/sqrt3; - q_Bragg[2]=ql_Bragg*cstar; + qsquare = 1; + Gn = 0; - for (i=0; i<4; i++) // Calculate phonon structure factor; loop over atoms in unit cell + double q_Bragg[3]; // Bragg point in Å-1; Cartesian coordinates + q_Bragg[0] = qh_Bragg * astar + qk_Bragg * astar / 2.0; + q_Bragg[1] = qk_Bragg * astar * 2.0 / sqrt3; + q_Bragg[2] = ql_Bragg * cstar; + + for (i = 0; i < 4; i++) // Calculate phonon structure factor; loop over atoms in unit cell + { + q_dot_e = (complex)0; + q_dot_r = 0; + for (j = 0; j < 3; j++) // loop over x,y,z { - q_dot_e = (complex)0; - q_dot_r = 0; - for (j=0; j<3; j++) // loop over x,y,z - { - q_dot_e += (complex)q[j]*eigenvectormode[3*i+j]; - q_dot_r += q_Bragg[j]*Delta[3*i+j]; - } - if (verbose >= 7) - { - printf("i = %i , q dot r = %g , q dot e = (%g + i %g) Fn contrib = (%g + i %g)",i,q_dot_r,creal(q_dot_e),cimag(q_dot_e),creal(b_length*q_dot_e*cexp(I*q_dot_r)),cimag(b_length*q_dot_e*cexp(I*q_dot_r))); - } - Gn += b_length*q_dot_e*cexp(I*q_dot_r); // Unit fm Å-1 ; for the general lattice, this should be divided by sqrt(mass[j]) + q_dot_e += (complex)q[j] * eigenvectormode[3 * i + j]; + q_dot_r += q_Bragg[j] * Delta[3 * i + j]; + } + if (verbose >= 7) { + printf ("i = %i , q dot r = %g , q dot e = (%g + i %g) Fn contrib = (%g + i %g)", i, q_dot_r, creal (q_dot_e), cimag (q_dot_e), + creal (b_length * q_dot_e * cexp (I * q_dot_r)), cimag (b_length * q_dot_e * cexp (I * q_dot_r))); } - if (verbose >= 7) - printf("\n"); + Gn += b_length * q_dot_e * cexp (I * q_dot_r); // Unit fm Å-1 ; for the general lattice, this should be divided by sqrt(mass[j]) + } + if (verbose >= 7) + printf ("\n"); if (shape == 0) - intersect = cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius, yheight); + intersect = cylinder_intersect (&t0, &t3, x, y, z, vx, vy, vz, radius, yheight); else if (shape == 1) - intersect = box_intersect(&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth); - - if(!intersect) - { - printf("FATAL ERROR: Did not hit sample surface from inside.\n"); - exit(1); - } - - if (verbose>=2) - printf("neutron left Phonon_BvK_PG, p_init = %g,",p); - - dt = t3; - l_o = v_f*dt; - my_a_i = V_my_a_v/v_i; - my_a_f = V_my_a_v/v_f; - bose_factor=nbose(omega,T); - J0 = hbar*hbar*v_f*V2K*1E10/mn; // The J-factor for a flat mode - - - p_att = exp(-(V_my_s*(l_i+l_o)+my_a_i*l_i+my_a_f*l_o)); /* Absorption factor (units: 1) */ - p_MC = nf*solid_angle*l_full; /* Focusing factors; assume random choice of n_f possibilities (units: m) */ - if (mode_input == 12) - p_MC *= 12; //The factor 12 is due to MC choice between the 12 modes - p_ph1 = (v_f/v_i)*DW*bose_factor*V_rho*1E30; /* Cross section factor 1 (units: m^-3) */ - p_ph2 = hbar*hbar/(2*(fabs(omega)*1.6022E-22)); /* Jacobian of delta functions in cross section; and constants. Units: J s^2 */ - p_ph3 = (Gn*conj(Gn))*1E-10/M; /* Structure factor. (units: kg^-1) */ - p_J = J0/J_factor; - p_tot = p_att*p_MC*p_ph1*p_ph2*p_ph3*p_J; - p *= p_tot; - - if (verbose>=2) { - printf("nf = %d, solid_angle = %g, l_full = %g \n",nf,solid_angle,l_full); - printf("(p in SI units): p_att = %g, p_MC= %g, p_ph1 = %g, p_ph2 = %g, p_ph3 = %g, p_J = %g, bose_factor = %g, V_rho = %g, J_factor/J0 = %g, G_factor = %g +i %g, conj(G) = %g + I %g, |G|^2 = %g + I %g, ",p_att, p_MC, p_ph1, p_ph2, p_ph3, p_J, bose_factor, V_rho, J_factor/J0, creal(Gn), cimag(Gn), creal(conj(Gn)), cimag(conj(Gn)), creal(Gn*conj(Gn)), cimag(Gn*conj(Gn)) ); - printf(" p_tot = %g, p_final = %g\n",p_tot, p) ; - printf("q is (h,k,l) = (%g, %g, %g) \n",q_x/astar,q_y/astar,q_z/cstar); - } + intersect = box_intersect (&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth); + + if (!intersect) { + printf ("FATAL ERROR: Did not hit sample surface from inside.\n"); + exit (1); + } + + if (verbose >= 2) + printf ("neutron left Phonon_BvK_PG, p_init = %g,", p); + + dt = t3; + l_o = v_f * dt; + my_a_i = V_my_a_v / v_i; + my_a_f = V_my_a_v / v_f; + bose_factor = nbose (omega, T); + J0 = hbar * hbar * v_f * V2K * 1E10 / mn; // The J-factor for a flat mode + + p_att = exp (-(V_my_s * (l_i + l_o) + my_a_i * l_i + my_a_f * l_o)); /* Absorption factor (units: 1) */ + p_MC = nf * solid_angle * l_full; /* Focusing factors; assume random choice of n_f possibilities (units: m) */ + if (mode_input == 12) + p_MC *= 12; // The factor 12 is due to MC choice between the 12 modes + p_ph1 = (v_f / v_i) * DW * bose_factor * V_rho * 1E30; /* Cross section factor 1 (units: m^-3) */ + p_ph2 = hbar * hbar / (2 * (fabs (omega) * 1.6022E-22)); /* Jacobian of delta functions in cross section; and constants. Units: J s^2 */ + p_ph3 = (Gn * conj (Gn)) * 1E-10 / M; /* Structure factor. (units: kg^-1) */ + p_J = J0 / J_factor; + p_tot = p_att * p_MC * p_ph1 * p_ph2 * p_ph3 * p_J; + p *= p_tot; + + if (verbose >= 2) { + printf ("nf = %d, solid_angle = %g, l_full = %g \n", nf, solid_angle, l_full); + printf ("(p in SI units): p_att = %g, p_MC= %g, p_ph1 = %g, p_ph2 = %g, p_ph3 = %g, p_J = %g, bose_factor = %g, V_rho = %g, J_factor/J0 = %g, G_factor = " + "%g +i %g, conj(G) = %g + I %g, |G|^2 = %g + I %g, ", + p_att, p_MC, p_ph1, p_ph2, p_ph3, p_J, bose_factor, V_rho, J_factor / J0, creal (Gn), cimag (Gn), creal (conj (Gn)), cimag (conj (Gn)), + creal (Gn * conj (Gn)), cimag (Gn * conj (Gn))); + printf (" p_tot = %g, p_final = %g\n", p_tot, p); + printf ("q is (h,k,l) = (%g, %g, %g) \n", q_x / astar, q_y / astar, q_z / cstar); + } } /* else transmit: Neutron did not hit the sample */ %} MCDISPLAY %{ - if (radius > 0 && yheight) { /* cylinder */ - cylinder(0,0,0,radius,yheight,0, 0, 1, 0); - } - 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); + if (radius > 0 && yheight) { /* cylinder */ + cylinder (0, 0, 0, radius, yheight, 0, 0, 1, 0); + } 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 From 9ef6b2efbc562b226f4f8cba206c3755217801d0 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 25 Feb 2026 21:53:50 +0100 Subject: [PATCH 14/15] AI-assisted rewrite of central functions to use mccode-complex-lib wrappers. (Work remains, some of the comments need to go back in etc...) --- mcstas-comps/contrib/Phonon_BvK_PG.comp | 818 ++++++++++-------------- 1 file changed, 335 insertions(+), 483 deletions(-) diff --git a/mcstas-comps/contrib/Phonon_BvK_PG.comp b/mcstas-comps/contrib/Phonon_BvK_PG.comp index 45c896122..6b969d411 100644 --- a/mcstas-comps/contrib/Phonon_BvK_PG.comp +++ b/mcstas-comps/contrib/Phonon_BvK_PG.comp @@ -366,13 +366,13 @@ SHARE } /* TODO: (James Avery) Replace by standard complex.h operation */ - complex + cdouble complexexp_qdotr (double* q, double* rj) { // double qrjtest double qrj = *(q) * *(rj) + *(q + 1) * *(rj + 1) + *(q + 2) * *(rj + 2); // printf(" q= (%g %g %g), r = (%g %g %g) qrj = %g \n",q[0],q[1],q[2],rj[0],rj[1],rj[2],qrj); // return (cexp(qrj)); // KL comment 030125 WHY does this not run? - return (cos (qrj) + I * sin (qrj)); + return cplx (cos (qrj), sin (qrj)); } void @@ -395,387 +395,281 @@ SHARE //---------------------------- START CENTRAL FUNCTION OMEGA-Q ---------------------------- double - omega_q (complex* parms, complex* eigenvectorgood) { - - // _______________ DETERMINE Q ____________________________ - + omega_q (cdouble* parms, cdouble* eigenvectorgood) { double vi, vf, vv_x, vv_y, vv_z, vi_x, vi_y, vi_z; double q, qx, qy, qz, Jq, hw_phonon, hw_neutron; double h, k, l, hk_inplane, hk_outofplane; int disp; - for (int ii = 0; ii < 8; ii++) - if (isnan (creal (parms[ii])) || isnan (cimag (parms[ii]))) { - fprintf (stderr, "omega_q: Array parms contains a nan at position %d.\n", ii); - fprintf (stderr, "parms = ["); - for (int i = 0; i < 8; i++) - fprintf (stderr, "%g + i%g%s", creal (parms[i]), cimag (parms[i]), i + 1 < 8 ? ", " : "]\n"); - abort (); - } + vf = creal (parms[0]); + vi = creal (parms[1]); + vv_x = creal (parms[2]); + vv_y = creal (parms[3]); + vv_z = creal (parms[4]); + vi_x = creal (parms[5]); + vi_y = creal (parms[6]); + vi_z = creal (parms[7]); - vf = parms[0]; - vi = parms[1]; - vv_x = parms[2]; - vv_y = parms[3]; - vv_z = parms[4]; - vi_x = parms[5]; - vi_y = parms[6]; - vi_z = parms[7]; - h = parms[9]; - k = parms[10]; - l = parms[11]; - disp = (int)parms[12]; + h = creal (parms[9]); + k = creal (parms[10]); + l = creal (parms[11]); + disp = (int)creal (parms[12]); qx = V2K * (vi_x - vf * vv_x); qy = V2K * (vi_y - vf * vv_y); qz = V2K * (vi_z - vf * vv_z); - if (disp == 1) /* calculate dispersion directly */ - { + if (disp == 1) { printf ("Dispersion calculation: "); qx = h * astar; - qy = k * astar; // correct only for qy=0 ! + qy = k * astar; qz = l * cstar; } double qvec[3] = { qx, qy, qz }; - if (verbose >= 4) - printf ("Enter omega_q; q = (%g, %g, %g) \n", qx, qy, qz); + /* ===================== Phi_diag ===================== */ - // ________________ END DETERMINE Q ________________- + cdouble Phi_diag[9]; - complex Phi_diag[3 * 3]; - int i; - int j; - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - Phi_diag[i * 3 + j] = Phi1[i][j] + Phi2[i][j] + Phi3[i][j] + Phi4[i][j] * (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[3][0])) - + Phi5[i][j] * (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[4][0])) + Phi6[i][j] * (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[5][0])) - + Phi7[i][j] * (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[6][0])) + Phi8[i][j] * (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[7][0])) - + Phi9[i][j] * (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[8][0])) + Phi10[i][j] + Phi11[i][j] + Phi12[i][j]; - } - } + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) { + cdouble sum = cplx (Phi1[i][j] + Phi2[i][j] + Phi3[i][j] + Phi10[i][j] + Phi11[i][j] + Phi12[i][j], 0.0); - complex Phi_diag1[3 * 3]; - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - Phi_diag1[i * 3 + j] = Phi13[i][j] + Phi14[i][j]; + cdouble e; + + e = (*complexexp_qdotr) (&qvec[0], &r_j13[3][0]); + sum = cadd (sum, rmul (Phi4[i][j], csub (cplx (1.0, 0.0), e))); + + e = (*complexexp_qdotr) (&qvec[0], &r_j13[4][0]); + sum = cadd (sum, rmul (Phi5[i][j], csub (cplx (1.0, 0.0), e))); + + e = (*complexexp_qdotr) (&qvec[0], &r_j13[5][0]); + sum = cadd (sum, rmul (Phi6[i][j], csub (cplx (1.0, 0.0), e))); + + e = (*complexexp_qdotr) (&qvec[0], &r_j13[6][0]); + sum = cadd (sum, rmul (Phi7[i][j], csub (cplx (1.0, 0.0), e))); + + e = (*complexexp_qdotr) (&qvec[0], &r_j13[7][0]); + sum = cadd (sum, rmul (Phi8[i][j], csub (cplx (1.0, 0.0), e))); + + e = (*complexexp_qdotr) (&qvec[0], &r_j13[8][0]); + sum = cadd (sum, rmul (Phi9[i][j], csub (cplx (1.0, 0.0), e))); + + Phi_diag[i * 3 + j] = sum; } - } - complex Phi_offdiag1[3 * 3]; - complex Phi_offdiag2[3 * 3]; - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - Phi_offdiag1[i * 3 + j] = -Phi1[i][j] * ((*complexexp_qdotr) (&qvec[0], &r_j13[0][0])) - Phi2[i][j] * ((*complexexp_qdotr) (&qvec[0], &r_j13[1][0])) - - Phi3[i][j] * ((*complexexp_qdotr) (&qvec[0], &r_j13[2][0])) - Phi10[i][j] * ((*complexexp_qdotr) (&qvec[0], &r_j13[9][0])) - - Phi11[i][j] * ((*complexexp_qdotr) (&qvec[0], &r_j13[10][0])) - Phi12[i][j] * ((*complexexp_qdotr) (&qvec[0], &r_j13[11][0])); - Phi_offdiag2[i * 3 + j] = conj (Phi_offdiag1[i * 3 + j]); + /* ===================== Phi_diag1 ===================== */ + + cdouble Phi_diag1[9]; + + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + Phi_diag1[i * 3 + j] = cplx (Phi13[i][j] + Phi14[i][j], 0.0); + + /* ===================== Phi_offdiag ===================== */ + + cdouble Phi_offdiag1[9]; + cdouble Phi_offdiag2[9]; + + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) { + cdouble s = cplx (0.0, 0.0); + + s = cadd (s, cneg (rmul (Phi1[i][j], (*complexexp_qdotr) (&qvec[0], &r_j13[0][0])))); + + s = cadd (s, cneg (rmul (Phi2[i][j], (*complexexp_qdotr) (&qvec[0], &r_j13[1][0])))); + + s = cadd (s, cneg (rmul (Phi3[i][j], (*complexexp_qdotr) (&qvec[0], &r_j13[2][0])))); + + s = cadd (s, cneg (rmul (Phi10[i][j], (*complexexp_qdotr) (&qvec[0], &r_j13[9][0])))); + + s = cadd (s, cneg (rmul (Phi11[i][j], (*complexexp_qdotr) (&qvec[0], &r_j13[10][0])))); + + s = cadd (s, cneg (rmul (Phi12[i][j], (*complexexp_qdotr) (&qvec[0], &r_j13[11][0])))); + + Phi_offdiag1[i * 3 + j] = s; + Phi_offdiag2[i * 3 + j] = conj (s); } - } - complex Phi_6Dim[6 * 6]; - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - Phi_6Dim[i * 6 + j] = Phi_diag[i * 3 + j] + Phi_diag1[i * 3 + j]; + /* ===================== Phi_6Dim ===================== */ + + cdouble Phi_6Dim[36]; + + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) { + Phi_6Dim[i * 6 + j] = cadd (Phi_diag[i * 3 + j], Phi_diag1[i * 3 + j]); Phi_6Dim[i * 6 + j + 3] = Phi_offdiag1[i * 3 + j]; Phi_6Dim[(i + 3) * 6 + j] = Phi_offdiag2[i * 3 + j]; - Phi_6Dim[(i + 3) * 6 + j + 3] = Phi_diag[i * 3 + j]; // Because Phi_diag1 does not appear for the B,C, sublattices + Phi_6Dim[(i + 3) * 6 + j + 3] = Phi_diag[i * 3 + j]; } - } - complex Phi_AC[3 * 3]; - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - Phi_AC[i * 3 + j] = -Phi13[i][j] * ((*complexexp_qdotr) (&qvec[0], &r_j13[12][0])) - Phi14[i][j] * ((*complexexp_qdotr) (&qvec[0], &r_j13[13][0])); + /* ===================== Phi_AC ===================== */ + + cdouble Phi_AC[9]; + + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) { + cdouble s = cplx (0.0, 0.0); + + s = cadd (s, cneg (rmul (Phi13[i][j], (*complexexp_qdotr) (&qvec[0], &r_j13[12][0])))); + + s = cadd (s, cneg (rmul (Phi14[i][j], (*complexexp_qdotr) (&qvec[0], &r_j13[13][0])))); + + Phi_AC[i * 3 + j] = s; } - } - // complex Zero_3Dim[3][3] = {{0 , 0 , 0},{0 , 0 , 0},{0 , 0 , 0}}; + /* ===================== Phi_6Dim_offdiag ===================== */ - complex Phi_6Dim_offdiag[6 * 6]; - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { + cdouble Phi_6Dim_offdiag[36]; + + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) { Phi_6Dim_offdiag[i * 6 + j] = Phi_AC[i * 3 + j]; - Phi_6Dim_offdiag[i * 6 + j + 3] = (complex)0.0; - Phi_6Dim_offdiag[(i + 3) * 6 + j] = (complex)0.0; - Phi_6Dim_offdiag[(i + 3) * 6 + j + 3] = (complex)0.0; + Phi_6Dim_offdiag[i * 6 + j + 3] = cplx (0.0, 0.0); + Phi_6Dim_offdiag[(i + 3) * 6 + j] = cplx (0.0, 0.0); + Phi_6Dim_offdiag[(i + 3) * 6 + j + 3] = cplx (0.0, 0.0); } - } - complex Phi_12Dim[12 * 12]; - for (i = 0; i < 6; i++) { - for (j = 0; j < 6; j++) { + /* ===================== Phi_12Dim ===================== */ + + cdouble Phi_12Dim[144]; + + for (int i = 0; i < 6; i++) + for (int j = 0; j < 6; j++) { Phi_12Dim[i * DIM + j] = Phi_6Dim[i * 6 + j]; Phi_12Dim[i * DIM + j + 6] = Phi_6Dim_offdiag[i * 6 + j]; Phi_12Dim[(i + 6) * DIM + j] = conj (Phi_6Dim_offdiag[i * 6 + j]); Phi_12Dim[(i + 6) * DIM + j + 6] = Phi_6Dim[i * 6 + j]; } - } - #ifdef TEST_MATRICES - printf ("Phi_12Dim=\n"); - for (i = 0; i < 12; i++) { - for (j = 0; j < 12; j++) { - if (j == 0) { - printf ("{"); - } - printf ("%g+(%g)i ", creal (Phi_12Dim[i * DIM + j]), cimag ((Phi_12Dim[i * DIM + j]))); - if (j == 11) { - printf ("}\n"); - } - } - } - - printf ("Phi_6Dim=\n"); - for (i = 0; i < 6; i++) { - for (j = 0; j < 6; j++) { - if (j == 0) { - printf ("{"); - } - printf ("%g+(%g)i ", creal (Phi_6Dim[i * 6 + j]), cimag ((Phi_6Dim[i * 6 + j]))); - if (j == 5) { - printf ("}\n"); - } - } - } - - printf ("Phi_diag=\n"); - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - if (j == 0) { - printf ("{"); - } - printf ("%g+(%g)i ", creal (Phi_diag[i * 3 + j]), cimag ((Phi_diag[i * 3 + j]))); - if (j == 2) { - printf ("}\n"); - } - } - } - printf ("Phi2={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}}\n", Phi2[0][0], Phi2[0][1], Phi2[0][2], Phi2[1][0], Phi2[1][1], Phi2[1][2], Phi2[2][0], Phi2[2][1], - Phi2[2][2]); - printf ("Phi3={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}", Phi3[0][0], Phi3[0][1], Phi3[0][2], Phi3[1][0], Phi3[1][1], Phi3[1][2], Phi3[2][0], Phi3[2][1], - Phi3[2][2]); - printf ("Phi5={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}", Phi5[0][0], Phi5[0][1], Phi5[0][2], Phi5[1][0], Phi5[1][1], Phi5[1][2], Phi5[2][0], Phi5[2][1], - Phi5[2][2]); - printf ("Phi6={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}", Phi6[0][0], Phi6[0][1], Phi6[0][2], Phi6[1][0], Phi6[1][1], Phi6[1][2], Phi6[2][0], Phi6[2][1], - Phi6[2][2]); - printf ("Phi7={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}", Phi7[0][0], Phi7[0][1], Phi7[0][2], Phi7[1][0], Phi7[1][1], Phi7[1][2], Phi7[2][0], Phi7[2][1], - Phi7[2][2]); - printf ("Phi8={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}", Phi8[0][0], Phi8[0][1], Phi8[0][2], Phi8[1][0], Phi8[1][1], Phi8[1][2], Phi8[2][0], Phi8[2][1], - Phi8[2][2]); - printf ("Phi9={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}", Phi9[0][0], Phi9[0][1], Phi9[0][2], Phi9[1][0], Phi9[1][1], Phi9[1][2], Phi9[2][0], Phi9[2][1], - Phi9[2][2]); - printf ("Phi11={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}", Phi11[0][0], Phi11[0][1], Phi11[0][2], Phi11[1][0], Phi11[1][1], Phi11[1][2], Phi11[2][0], - Phi11[2][1], Phi11[2][2]); - printf ("Phi12={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}", Phi12[0][0], Phi12[0][1], Phi12[0][2], Phi12[1][0], Phi12[1][1], Phi12[1][2], Phi12[2][0], - Phi12[2][1], Phi12[2][2]); - printf ("Phi14={{%lf,%lf,%lf},{%lf,%lf,%lf},{%lf,%lf,%lf}\n}", Phi14[0][0], Phi14[0][1], Phi14[0][2], Phi14[1][0], Phi14[1][1], Phi14[1][2], Phi14[2][0], - Phi14[2][1], Phi14[2][2]); - - printf (" 1-exp(q.r_j13[3]) = %lf + i %lf \n", creal (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[3][0])), - cimag (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[3][0]))); - printf (" 1-exp(q.r_j13[4]) = %lf + i %lf \n", creal (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[4][0])), - cimag (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[4][0]))); - printf (" 1-exp(q.r_j13[5]) = %lf + i %lf \n", creal (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[5][0])), - cimag (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[5][0]))); - printf (" 1-exp(q.r_j13[6]) = %lf + i %lf \n", creal (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[6][0])), - cimag (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[6][0]))); - printf (" 1-exp(q.r_j13[7]) = %lf + i %lf \n", creal (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[7][0])), - cimag (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[7][0]))); - printf (" 1-exp(q.r_j13[8]) = %lf + i %lf \n", creal (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[8][0])), - cimag (1 - (*complexexp_qdotr) (&qvec[0], &r_j13[8][0]))); - - #endif // TEST_MATRICES + /* ===================== Build Hermitian matrix ===================== */ double eigenvalue[DIM]; - complex Q[DIM * DIM]; - complex matrix[DIM * DIM]; + cdouble Q[DIM * DIM]; + cdouble matrix[DIM * DIM]; int index[DIM]; - for (i = 0; i < DIM; i++) // transforms from 2D matrix to 1D pseudo-matrix - for (j = i; j < DIM; j++) { + for (int i = 0; i < DIM; i++) + for (int j = i; j < DIM; j++) { matrix[i * DIM + j] = Phi_12Dim[i * DIM + j]; - matrix[j * DIM + i] = conj (matrix[i * DIM + j]); // Forces the pseudo-matrix Hermitean - } // m(12*i+j) = M(i,j) - - #if MATRIX_TEST - fwrite (matrix, sizeof (complex), DIM * DIM, matrix_test_file); - fwrite (parms, sizeof (complex), 8, parms_test_file); - #endif + matrix[j * DIM + i] = conj (matrix[i * DIM + j]); + } - eigensystem_hermitian (matrix, DIM, eigenvalue, Q); /* Produces unitary transformation matrix Q: eigenvectors are columns */ + eigensystem_hermitian (matrix, DIM, eigenvalue, Q); - for (i = 0; i < DIM; i++) + for (int i = 0; i < DIM; i++) index[i] = i; - if (verbose >= 6) { - printf ("Before bubblesort\n"); - for (i = 0; i < DIM; i++) { - printf ("i = %i eigenvalue[i] = %g index[i] %i \n", i, eigenvalue[i], index[i]); - } - } - - // Sort eigenvalues from lowest to highest - bubblesort (DIM, eigenvalue, index); // mergesort(eigenvalue,DIM,index); - if (verbose >= 5) { - printf ("After bubblesort\n"); - for (i = 0; i < DIM; i++) { - printf ("i = %i eigenvalue[i] = %g index[i] %i \n", i, eigenvalue[i], index[i]); - } - } + bubblesort (DIM, eigenvalue, index); double eigenvaluegood[DIM]; - for (j = 0; j < DIM; j++) { - eigenvectorgood[j] = Q[index[mode] * DIM + j]; /* Eigenvectors are columns of Q */ // No, they are rows! - } + for (int j = 0; j < DIM; j++) + eigenvectorgood[j] = Q[index[mode] * DIM + j]; - if (verbose >= 7) { - printf ("Q = ["); - for (i = 0; i < DIM; i++) { - printf ("\n ("); - for (j = 0; j < DIM; j++) - printf (" %g +i %g, ", creal (Q[i * DIM + j]), cimag (Q[i * DIM + j])); - printf (")"); - } - printf ("]\n"); - } + for (int i = 0; i < DIM; i++) + eigenvaluegood[i] = sqrt (eigenvalue[index[i]]) / sqrt (M) / (2 * PI * 1E12) * THz2meV; - for (i = 0; i < DIM; i++) - eigenvaluegood[i] = sqrt (creal (eigenvalue[index[i]])) / sqrt (M) / (2 * PI * 1E12) * THz2meV; // convert to units of THz - - if (verbose >= 5) { - printf ("Diagonalization done, eigenvalues (energies) are: \n ("); - for (i = 0; i < DIM; i++) - printf (" %g, ", eigenvaluegood[i]); - printf (" )\n"); - printf ("mode = %i , eigenvectorgood[mode] = (", mode); - for (j = 0; j < DIM; j++) - printf (" %g +i %g, ", creal (eigenvectorgood[j]), cimag (eigenvectorgood[j])); - printf (" )\n"); - } + hw_neutron = fabs (VS2E * (vi * vi - vf * vf)); + hw_phonon = eigenvaluegood[mode]; - // ___________________ RETURN RESULTS TO CALLING FUNCTION __________________ - - hw_neutron = fabs (VS2E * (vi * vi - vf * vf)); // neutron energy transfer - // fabs used to find both positive and negative solutions, controlled by findroots() - hw_phonon = eigenvaluegood[mode]; // phonon energy - - if (verbose >= 4) { - printf ("Return from omega_q \n"); - printf ("hw_neutron = %g (vi = %g, vf = %g)\n", hw_neutron, vi, vf); - printf ("hw_phonon = %g = eigenvaluegood[%d] = %g = eigenvalues[%d] = %g\n", hw_phonon, mode, eigenvaluegood[mode], index[mode], eigenvalue[index[mode]]); - // printf("old parms = ["); for(int i=0;i<8;i++) fprintf(stderr,"%g+i%g%s",creal(parms[i]), cimag(parms[i]), i+1<8?", ":"]\n"); - } - - parms[8] = hw_phonon; - - if (disp == 1) { - printf ("mode = %d, (h,k,l) = (%g,%g,%g), hw_q = %g", mode, h, k, l, hw_phonon); - printf (" polarization: ( "); - for (j = 0; j < DIM; j++) { - printf ("%g + i %g ,", creal (eigenvectorgood[j]), cimag (eigenvectorgood[j])); - } - printf (")\n"); - } - if (disp == 2) { - hk_inplane = qx / astar; - hk_outofplane = qy / astar; // Derivation assume astar along x and 60 degrees between astar and bstar: - // h_ip = h + k/2 h_op = sqrt(3)*k/2 => k = 2/sqrt(3)*h_op h = h_ip - h_op/sqrt(3) - k = 2 / sqrt (3) * hk_outofplane; - h = hk_inplane - hk_outofplane / sqrt (3); - l = qz / cstar; - printf ("mode = %d, vf = %g, (h,k,l) = (%g,%g,%g), hw_q = %g", mode, vf, h, k, l, hw_phonon); - printf (" polarization: (%g + i %g , %g + i %g , %g + i %g) \n", creal (eigenvectorgood[0]), cimag (eigenvectorgood[0]), creal (eigenvectorgood[1]), - cimag (eigenvectorgood[1]), creal (eigenvectorgood[2]), cimag (eigenvectorgood[2])); - } + parms[8] = cplx (hw_phonon, 0.0); return (hw_phonon - hw_neutron); } - //--------------------- END CENTRAL FUNCTION OMEGA-Q ------------------------ // ------------------START OF THE GENERAL ROOT-FINDING CODE-------------------- double last_fh, last_x2; double - zridd (double (*func) (complex*, complex*), double x1, double x2, complex* parms, complex* vector, double xacc) { + zridd (double (*func) (cdouble*, cdouble*), double x1, double x2, cdouble* parms, cdouble* vector, double xacc) { int j; double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; - // fprintf(stderr,"i zridd(%g,%g, parms,%g)\n",x1,x2,xacc); - parms[0] = x1; - // fprintf(stderr,"zridd fl = omega_q\n"); + /* ---- f(x1) ---- */ + + parms[0] = cplx (x1, 0.0); + if (xacc > 0) { fl = (*func) (parms, vector); } else { fl = last_fh; xacc = -xacc; - if (fabs (x1 - last_x2) > xacc) // KL comment 030125 WHY does the library routine zridd use complex abs() function ?? changed to fabs() + + if (fabs (x1 - last_x2) > xacc) printf ("*** error in zridd *** x1: %g last_x2: %g \n", x1, last_x2); } - // fprintf(stderr,"fl = %g\n",fl); - parms[0] = x2; - // fprintf(stderr,"zridd fh = omega_q; parms[0] = x2 = %g+i%g = %g\n",creal(parms[0]), cimag(parms[0]),x2); + + /* ---- f(x2) ---- */ + + parms[0] = cplx (x2, 0.0); last_fh = fh = (*func) (parms, vector); last_x2 = x2; - // fprintf(stderr,"fh = %g\n",fh); + if (fl * fh >= 0) { if (fl == 0) return x1; if (fh == 0) return x2; return UNUSED; - } else { - xl = x1; - xh = x2; - // printf("zridd sign change: v_low : %g v_high: %g \n",xl,xh); - ans = UNUSED; - for (j = 1; j < MAXRIDD; j++) { - xm = 0.5 * (xl + xh); - parms[0] = xm; - // fprintf(stderr,"zridd fm = omega_q\n"); - fm = (*func) (parms, vector); - s = sqrt (fm * fm - fl * fh); - if (s == 0.0) - return ans; - xnew = xm + (xm - xl) * ((fl >= fh ? 1.0 : -1.0) * fm / s); - if (fabs (xnew - ans) <= xacc) - return ans; - ans = xnew; - parms[0] = ans; - // fprintf(stderr,"s = %g; fl, fm, fh = %g,%g,%g; fnew = %g; fm*fm-fl*fh = %g\n", - // s, fl, fm,fh, fnew, fm*fm-fl*fh); - // fprintf(stderr,"zridd fnew = omega_q\n"); - fnew = (*func) (parms, vector); - if (fnew == 0.0) - return ans; - if (fabs (fm) * SIGN (fnew) != fm) { - xl = xm; - fl = fm; - xh = ans; - fh = fnew; - } else if (fabs (fl) * SIGN (fnew) != fl) { - xh = ans; - fh = fnew; - } else if (fabs (fh) * SIGN (fnew) != fh) { - xl = ans; - fl = fnew; - } else - fatalerror ("never get here in zridd"); - if (fabs (xh - xl) <= xacc) - return ans; + } + + xl = x1; + xh = x2; + ans = UNUSED; + + for (j = 1; j < MAXRIDD; j++) { + xm = 0.5 * (xl + xh); + + parms[0] = cplx (xm, 0.0); + fm = (*func) (parms, vector); + + s = sqrt (fm * fm - fl * fh); + if (s == 0.0) + return ans; + + xnew = xm + (xm - xl) * ((fl >= fh ? 1.0 : -1.0) * fm / s); + + if (fabs (xnew - ans) <= xacc) + return ans; + + ans = xnew; + + parms[0] = cplx (ans, 0.0); + fnew = (*func) (parms, vector); + + if (fnew == 0.0) + return ans; + + if (fabs (fm) * SIGN (fnew) != fm) { + xl = xm; + fl = fm; + xh = ans; + fh = fnew; + } else if (fabs (fl) * SIGN (fnew) != fl) { + xh = ans; + fh = fnew; + } else if (fabs (fh) * SIGN (fnew) != fh) { + xl = ans; + fl = fnew; + } else { + fatalerror ("never get here in zridd"); } - fatalerror ("zridd exceeded maximum iterations"); + + if (fabs (xh - xl) <= xacc) + return ans; } - return 0.0; /* Never get here */ + + fatalerror ("zridd exceeded maximum iterations"); + return 0.0; } + #ifdef OPENACC #pragma acc routine double - zridd_gpu (double x1, double x2, complex* parms, complex* vector, double xacc) { + zridd_gpu (double x1, double x2, cdouble* parms, cdouble* vector, double xacc) { int j; double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew; @@ -832,29 +726,29 @@ SHARE } return 0.0; /* Never get here */ } + #endif #define ROOTACC 1e-8 int - findroots (double brack_low, double brack_mid, double brack_high, double* list, int* index, double (*f) (complex*, complex*), complex* parms, complex* vector, - int low_steps, int high_steps) - // *f is here omega_q - // *parms is here p_call - { + findroots (double brack_low, double brack_mid, double brack_high, double* list, int* index, double (*f) (cdouble*, cdouble*), cdouble* parms, cdouble* vector, + int low_steps, int high_steps) { double root, range; int i; if (verbose == 2) printf ("Enter findroots() \n"); - // First, find roots in energy loss side, if any + /* ---- Energy loss side ---- */ + range = brack_mid - brack_low; + for (i = 0; i < (low_steps - 1); i++) { if (i == 0) - root = zridd (f, brack_low + range * i / low_steps, brack_low + range * (i + 1) / low_steps, (complex*)parms, (complex*)vector, ROOTACC); + root = zridd (f, brack_low + range * i / low_steps, brack_low + range * (i + 1) / low_steps, parms, vector, ROOTACC); else - root = zridd (f, brack_low + range * i / low_steps, brack_low + range * (i + 1) / low_steps, (complex*)parms, (complex*)vector, - -ROOTACC); // Re-use the previous fh value + root = zridd (f, brack_low + range * i / low_steps, brack_low + range * (i + 1) / low_steps, parms, vector, -ROOTACC); + if (root != UNUSED) { list[(*index)++] = root; if (verbose >= 4) @@ -863,14 +757,15 @@ SHARE } range = (brack_mid - brack_low) / (double)low_steps; - for (i = 0; i < low_steps; i++) // Small steps close to zero energy transfer - { - if (low_steps == 1) // then the previous loop did not execute - root = zridd (f, brack_low + range * (low_steps - 1 + i / (double)low_steps), brack_low + range * (low_steps - 1 + (i + 1) / (double)low_steps), - (complex*)parms, (complex*)vector, ROOTACC); + + for (i = 0; i < low_steps; i++) { + if (low_steps == 1) + root = zridd (f, brack_low + range * (low_steps - 1 + i / (double)low_steps), brack_low + range * (low_steps - 1 + (i + 1) / (double)low_steps), parms, + vector, ROOTACC); else - root = zridd (f, brack_low + range * (low_steps - 1 + i / (double)low_steps), brack_low + range * (low_steps - 1 + (i + 1) / (double)low_steps), - (complex*)parms, (complex*)vector, -ROOTACC); // Re-use the previous fh value + root = zridd (f, brack_low + range * (low_steps - 1 + i / (double)low_steps), brack_low + range * (low_steps - 1 + (i + 1) / (double)low_steps), parms, + vector, -ROOTACC); + if (root != UNUSED) { list[(*index)++] = root; if (verbose >= 4) @@ -878,16 +773,16 @@ SHARE } } - // Second, find roots in energy gain side, there is always some + /* ---- Energy gain side ---- */ + range = (brack_high - brack_mid) / (double)high_steps; - for (i = 0; i < high_steps; i++) // Close to zero: small steps - { + + for (i = 0; i < high_steps; i++) { if (i == 0) - root - = zridd (f, brack_mid + range * i / (double)high_steps, brack_mid + range * (i + 1) / (double)high_steps, (complex*)parms, (complex*)vector, ROOTACC); + root = zridd (f, brack_mid + range * i / (double)high_steps, brack_mid + range * (i + 1) / (double)high_steps, parms, vector, ROOTACC); else - root = zridd (f, brack_mid + range * i / (double)high_steps, brack_mid + range * (i + 1) / (double)high_steps, (complex*)parms, (complex*)vector, - -ROOTACC); // Re-use the previous fh value + root = zridd (f, brack_mid + range * i / (double)high_steps, brack_mid + range * (i + 1) / (double)high_steps, parms, vector, -ROOTACC); + if (root != UNUSED) { list[(*index)++] = root; if (verbose >= 4) @@ -896,37 +791,41 @@ SHARE } range = brack_high - brack_mid; - for (i = 1; i < high_steps; i++) // Larger steps away from zero - { - root = zridd (f, brack_mid + range * i / (double)high_steps, brack_mid + range * (i + 1) / (double)high_steps, (complex*)parms, (complex*)vector, - -ROOTACC); // Re-use the previous fh value - there was always one + + for (i = 1; i < high_steps; i++) { + root = zridd (f, brack_mid + range * i / (double)high_steps, brack_mid + range * (i + 1) / (double)high_steps, parms, vector, -ROOTACC); + if (root != UNUSED) { list[(*index)++] = root; if (verbose >= 4) printf ("*** findroots returned a root on the energy gain side: vf = %g \n", root); } } + + return *index; } + #ifdef OPENACC #pragma acc routine int - findroots_gpu (double brack_low, double brack_mid, double brack_high, double* list, int* index, double* parms, complex* vector, int low_steps, int high_steps) { + findroots_gpu (double brack_low, double brack_mid, double brack_high, double* list, int* index, double* parms, cdouble* vector, int low_steps, int high_steps) { double root, range = brack_mid - brack_low; int i; for (i = 0; i < low_steps; i++) { - root = zridd_gpu (brack_low + range * i / (int)low_steps, brack_low + range * (i + 1) / (int)low_steps, (complex*)parms, (complex*)vector, ROOTACC); + root = zridd_gpu (brack_low + range * i / (int)low_steps, brack_low + range * (i + 1) / (int)low_steps, (cdouble*)parms, (cdouble*)vector, ROOTACC); if (root != UNUSED) { list[(*index)++] = root; } } for (i = 0; i < high_steps; i++) { - root = zridd_gpu (brack_mid + range * i / (int)high_steps, brack_high + range * (i + 1) / (int)high_steps, (complex*)parms, (complex*)vector, ROOTACC); + root = zridd_gpu (brack_mid + range * i / (int)high_steps, brack_high + range * (i + 1) / (int)high_steps, (cdouble*)parms, (cdouble*)vector, ROOTACC); if (root != UNUSED) { list[(*index)++] = root; } } } + #endif #undef UNUSED #undef MAXRIDD @@ -938,7 +837,7 @@ DECLARE double V_rho; double V_my_s; double V_my_a_v; - complex** Matrix; + cdouble** Matrix; int i, j; double q[3]; /* Scattering vector */ double qx; @@ -947,8 +846,8 @@ DECLARE double q_x; double q_y; double q_z; - complex eigenvectormode[DIM]; - complex p_call[15]; /* Parameter list to transfer to omega_q; elsewhere known as parms */ + cdouble eigenvectormode[DIM]; + cdouble p_call[15]; /* Parameter list to transfer to omega_q; elsewhere known as parms */ %} INITIALIZE @@ -990,29 +889,29 @@ INITIALIZE TRACE %{ - double t0, t1, t2, t3; /* Entry/exit times for shape */ - double v_i, v_f; /* Neutron velocities: initial, final */ - double vx_i, vy_i, vz_i; /* Neutron initial velocity vector */ - double dt0, dt; /* Flight times through sample */ - double l_full; /* Flight path length for non-scattered neutron */ - double l_i, l_o; /* Flight path lenght in/out for scattered neutron */ - double my_a_i; /* Initial attenuation factor */ - double my_a_f; /* Final attenuation factor */ - double solid_angle; /* Solid angle of target as seen from scattering point */ - double aim_x = 0, aim_y = 0, aim_z = 1; /* Position of target relative to scattering point */ - double omega; /* Energy transfer */ - double qsquare; /* Square of the scattering vector */ - double bose_factor; /* Calculated value of the Bose factor */ - int nf, index; /* Number of allowed final velocities */ - double vf_list[20]; /* List of allowed final velocities */ - double J_factor, J0; /* Jacobian from delta fnc.s in cross section; elastic J-factor */ - double f1, f2; /* Probed values of omega_q minus omega */ - double p_att, p_MC, p_ph1, p_ph2, p_ph3, p_J, p_tot; /* Temporary multipliers */ - complex Gn; /* Phonon structure factor */ - complex q_dot_e; /* q dot e */ - double q_dot_r; /* q dot r */ - double qh, qk, ql, qhkx, qhky, qh_Bragg, qk_Bragg, ql_Bragg, qh_res, qk_res, qh_add, qk_add; /* components of q */ - double dist_00, dist_01, dist_10, dist_11, dist_min; /* used to calculate nearest Bragg point */ + double t0, t1, t2, t3; + double v_i, v_f; + double vx_i, vy_i, vz_i; + double dt0, dt; + double l_full; + double l_i, l_o; + double my_a_i; + double my_a_f; + double solid_angle; + double aim_x = 0, aim_y = 0, aim_z = 1; + double omega; + double qsquare; + double bose_factor; + int nf, index; + double vf_list[20]; + double J_factor, J0; + double f1, f2; + double p_att, p_MC, p_ph1, p_ph2, p_ph3, p_J, p_tot; + cdouble Gn; + cdouble q_dot_e; + double q_dot_r; + double qh, qk, ql, qhkx, qhky, qh_Bragg, qk_Bragg, ql_Bragg, qh_res, qk_res, qh_add, qk_add; + double dist_00, dist_01, dist_10, dist_11, dist_min; int i, j, intersect; if (shape == 0) @@ -1021,119 +920,89 @@ TRACE intersect = box_intersect (&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth); if (intersect) { - // fprintf(stderr,"intersect\n"); if (t0 < 0) - ABSORB; /* Neutron came from the sample or begins inside */ + ABSORB; if (verbose >= 2) - printf ("neutron entered Phonon_BvK_PG \n"); - /* Neutron enters at t=t0. */ - dt0 = t3 - t0; /* Time in sample */ + printf ("neutron entered Phonon_BvK_PG\n"); + + dt0 = t3 - t0; v_i = sqrt (vx * vx + vy * vy + vz * vz); - l_full = v_i * dt0; /* Length of path through sample if not scattered */ - dt = rand01 () * dt0; /* Time of scattering (relative to t0) */ - l_i = v_i * dt; /* Penetration in sample at scattering */ + l_full = v_i * dt0; + dt = rand01 () * dt0; + l_i = v_i * dt; + vx_i = vx; vy_i = vy; vz_i = vz; - PROP_DT (dt + t0); /* Point of scattering */ - aim_x = target_x - x; /* Vector pointing at target (e.g. analyzer) */ + PROP_DT (dt + t0); + + aim_x = target_x - x; aim_y = target_y - y; aim_z = target_z - z; - if (focus_aw && focus_ah) { + if (focus_aw && focus_ah) randvec_target_rect_angular (&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_aw, focus_ah, ROT_A_CURRENT_COMP); - } else if (focus_xw && focus_yh) { + else if (focus_xw && focus_yh) randvec_target_rect (&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP); - } else { + else randvec_target_sphere (&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_r); - } + NORM (vx, vy, vz); + nf = 0; + if (mode_input == 12) - mode = round (12 * rand01 () - 0.5); // Select the phonon mode + mode = round (12 * rand01 () - 0.5); else - mode = mode_input; // Use the given mode - if ((mode < 0) || (mode > 11)) // Undefined mode specification - { + mode = mode_input; + + if ((mode < 0) || (mode > 11)) { printf ("mode = %d ", mode); nrerror ("illegal value of mode"); } - p_call[0] = -1; // in the end this will be v_f - p_call[1] = v_i; - p_call[2] = vx; - p_call[3] = vy; - p_call[4] = vz; - p_call[5] = vx_i; - p_call[6] = vy_i; - p_call[7] = vz_i; - p_call[9] = hh; - p_call[10] = kk; - p_call[11] = ll; - p_call[12] = (double)dispersion; + p_call[0] = cplx (-1, 0); + p_call[1] = cplx (v_i, 0); + p_call[2] = cplx (vx, 0); + p_call[3] = cplx (vy, 0); + p_call[4] = cplx (vz, 0); + p_call[5] = cplx (vx_i, 0); + p_call[6] = cplx (vy_i, 0); + p_call[7] = cplx (vz_i, 0); + p_call[9] = cplx (hh, 0); + p_call[10] = cplx (kk, 0); + p_call[11] = cplx (ll, 0); + p_call[12] = cplx ((double)dispersion, 0); if (dispersion == 1) { f1 = omega_q (p_call, eigenvectormode); - p_call[12] = 0.0; + p_call[12] = cplx (0, 0); } - /* if (dispersion==2) - { - for(i=0; i= 2) - printf ("Call findroots \n"); findroots (0, v_i, v_i + V_HIGH, vf_list, &nf, omega_q, p_call, eigenvectormode, e_steps_low, e_steps_high); - if (verbose >= 2) - printf ("Findroots returned %d roots \n", nf); #else findroots_gpu (0, v_i, v_i + V_HIGH, vf_list, &nf, p_call, eigenvectormode, e_steps_low, e_steps_high); #endif - if (verbose >= 3) { - printf ("Findroots done, mode %i , last phonon energy %g \n", mode, creal (p_call[8])); - } - - // ________________ ALL FROM HERE IS CALCULATION OF INTENSITY - SHOULD BE REVISITED ___________________ - - index = (int)floor (rand01 () * nf); // Select random solution + index = (int)floor (rand01 () * nf); v_f = vf_list[index]; - p_call[0] = v_f; // transfer choice of v_f to omega_q + p_call[0] = cplx (v_f, 0); f1 = omega_q (p_call, eigenvectormode); - if (verbose >= 2) { - printf ("Chosen solution is %i; v_f = %g, hw = %g, energy match is %g; eigenvector is: (", index, v_f, creal (p_call[8]), f1); - for (j = 0; j < DIM; j++) - printf (" %g +i %g, ", creal (eigenvectormode[j]), cimag (eigenvectormode[j])); - printf (" )\n"); - printf ("--------------------------------------------\n"); - } - p_call[0] = v_f - DV; - // fprintf(stderr,"f1=omega_q(%g+i%g)\n",creal(p_call[0]),cimag(p_call[0])); + + p_call[0] = cplx (v_f - DV, 0); f1 = omega_q (p_call, eigenvectormode); - // fprintf(stderr,"retur fra f1=omega_q\n"); - p_call[0] = v_f + DV; - // fprintf(stderr,"f2=omega_q(%g+i%g)\n",creal(p_call[0]),cimag(p_call[0])); + + p_call[0] = cplx (v_f + DV, 0); f2 = omega_q (p_call, eigenvectormode); - // fprintf(stderr,"retur fra f2=omega_q\n"); - J_factor = 1.6022E-22 * 1E-10 * fabs (f2 - f1) / (2 * DV * V2K); // unit: meV Å converted to SI units - omega = VS2E * (v_i * v_i - v_f * v_f); // unit: meV + + J_factor = 1.6022E-22 * 1E-10 * fabs (f2 - f1) / (2 * DV * V2K); + + omega = VS2E * (v_i * v_i - v_f * v_f); + vx *= v_f; vy *= v_f; vz *= v_f; @@ -1143,80 +1012,67 @@ TRACE q_z = q[2] = V2K * (vz_i - vz); ql = q_z / cstar; - qhkx = q_x / astar; // These are the Cartesian coordinates - qhky = q_y / astar; // These are the Cartesian coordinates - qk = 2 * qhky / sqrt3; // transform to hexagonal coordinates - qh = qhkx - qhky / sqrt3; // transform to hexagonal coordinates - - // Now find out which Bragg peak is closest in hexagonal plane + qhkx = q_x / astar; + qhky = q_y / astar; + qk = 2 * qhky / sqrt3; + qh = qhkx - qhky / sqrt3; - ql_Bragg = round (ql); // We know what to do along the c-axis + ql_Bragg = round (ql); - qh_Bragg = floor (qh); // Divide into integer and residual part + qh_Bragg = floor (qh); qh_res = qh - qh_Bragg; qk_Bragg = floor (qk); qk_res = qk - qk_Bragg; - dist_00 = qh_res * qh_res + qk_res * qk_res + qh_res * qk_res; // Square distance in hexagonal coordinates to (0,0) + dist_00 = qh_res * qh_res + qk_res * qk_res + qh_res * qk_res; dist_min = dist_00; qh_add = 0; qk_add = 0; - dist_01 = qh_res * qh_res + (qk_res - 1) * (qk_res - 1) + qh_res * (qk_res - 1); // Square distance in hexagonal coordinates to (0,1) + dist_01 = qh_res * qh_res + (qk_res - 1) * (qk_res - 1) + qh_res * (qk_res - 1); if (dist_01 < dist_min) { dist_min = dist_01; qh_add = 0; qk_add = 1; } - dist_10 = (qh_res - 1) * (qh_res - 1) + qk_res * qk_res + (qh_res - 1) * qk_res; // Square distance in hexagonal coordinates to (1,0) + + dist_10 = (qh_res - 1) * (qh_res - 1) + qk_res * qk_res + (qh_res - 1) * qk_res; if (dist_10 < dist_min) { dist_min = dist_10; qh_add = 1; qk_add = 0; } - dist_11 = (qh_res - 1) * (qh_res - 1) + (qk_res - 1) * (qk_res - 1) + (qh_res - 1) * (qk_res - 1); // Square distance in hexagonal coordinates to (1,1) + + dist_11 = (qh_res - 1) * (qh_res - 1) + (qk_res - 1) * (qk_res - 1) + (qh_res - 1) * (qk_res - 1); if (dist_11 < dist_min) { dist_min = dist_11; qh_add = 1; qk_add = 1; } - qh_Bragg += qh_add; // Potentially correct to the right BZ + qh_Bragg += qh_add; qk_Bragg += qk_add; - if (verbose >= 3) { - printf ("q transforms in r.l.u to (h, k, l) = (%g, %g, %g). Nearest Bragg point: (%g, %g, %g)\n", qh, qk, ql, qh_Bragg, qk_Bragg, ql_Bragg); - printf ("--------------------------------------------\n"); - } - - if (dispersion == 1) - printf ("Dispersion found: (h,k,l) = (%g, %g, %g), hw = %g \n", q[0] / astar, q[1] / astar, q[2] / cstar, omega); - qsquare = 1; - Gn = 0; + Gn = cplx (0, 0); - double q_Bragg[3]; // Bragg point in Å-1; Cartesian coordinates + double q_Bragg[3]; q_Bragg[0] = qh_Bragg * astar + qk_Bragg * astar / 2.0; q_Bragg[1] = qk_Bragg * astar * 2.0 / sqrt3; q_Bragg[2] = ql_Bragg * cstar; - for (i = 0; i < 4; i++) // Calculate phonon structure factor; loop over atoms in unit cell - { - q_dot_e = (complex)0; + for (i = 0; i < 4; i++) { + q_dot_e = cplx (0, 0); q_dot_r = 0; - for (j = 0; j < 3; j++) // loop over x,y,z - { - q_dot_e += (complex)q[j] * eigenvectormode[3 * i + j]; + + for (j = 0; j < 3; j++) { + q_dot_e = cadd (q_dot_e, rmul (q[j], eigenvectormode[3 * i + j])); q_dot_r += q_Bragg[j] * Delta[3 * i + j]; } - if (verbose >= 7) { - printf ("i = %i , q dot r = %g , q dot e = (%g + i %g) Fn contrib = (%g + i %g)", i, q_dot_r, creal (q_dot_e), cimag (q_dot_e), - creal (b_length * q_dot_e * cexp (I * q_dot_r)), cimag (b_length * q_dot_e * cexp (I * q_dot_r))); - } - Gn += b_length * q_dot_e * cexp (I * q_dot_r); // Unit fm Å-1 ; for the general lattice, this should be divided by sqrt(mass[j]) + + cdouble phase = cexp (cplx (0, q_dot_r)); + Gn = cadd (Gn, cmul (cplx (b_length, 0), cmul (q_dot_e, phase))); } - if (verbose >= 7) - printf ("\n"); if (shape == 0) intersect = cylinder_intersect (&t0, &t3, x, y, z, vx, vy, vz, radius, yheight); @@ -1224,41 +1080,37 @@ TRACE intersect = box_intersect (&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth); if (!intersect) { - printf ("FATAL ERROR: Did not hit sample surface from inside.\n"); + printf ("FATAL ERROR\n"); exit (1); } - if (verbose >= 2) - printf ("neutron left Phonon_BvK_PG, p_init = %g,", p); - dt = t3; l_o = v_f * dt; + my_a_i = V_my_a_v / v_i; my_a_f = V_my_a_v / v_f; + bose_factor = nbose (omega, T); - J0 = hbar * hbar * v_f * V2K * 1E10 / mn; // The J-factor for a flat mode - p_att = exp (-(V_my_s * (l_i + l_o) + my_a_i * l_i + my_a_f * l_o)); /* Absorption factor (units: 1) */ - p_MC = nf * solid_angle * l_full; /* Focusing factors; assume random choice of n_f possibilities (units: m) */ + J0 = hbar * hbar * v_f * V2K * 1E10 / mn; + + p_att = exp (-(V_my_s * (l_i + l_o) + my_a_i * l_i + my_a_f * l_o)); + p_MC = nf * solid_angle * l_full; + if (mode_input == 12) - p_MC *= 12; // The factor 12 is due to MC choice between the 12 modes - p_ph1 = (v_f / v_i) * DW * bose_factor * V_rho * 1E30; /* Cross section factor 1 (units: m^-3) */ - p_ph2 = hbar * hbar / (2 * (fabs (omega) * 1.6022E-22)); /* Jacobian of delta functions in cross section; and constants. Units: J s^2 */ - p_ph3 = (Gn * conj (Gn)) * 1E-10 / M; /* Structure factor. (units: kg^-1) */ + p_MC *= 12; + + p_ph1 = (v_f / v_i) * DW * bose_factor * V_rho * 1E30; + p_ph2 = hbar * hbar / (2 * (fabs (omega) * 1.6022E-22)); + + cdouble G2 = cmul (Gn, conj (Gn)); + p_ph3 = creal (G2) * 1E-10 / M; + p_J = J0 / J_factor; + p_tot = p_att * p_MC * p_ph1 * p_ph2 * p_ph3 * p_J; p *= p_tot; - - if (verbose >= 2) { - printf ("nf = %d, solid_angle = %g, l_full = %g \n", nf, solid_angle, l_full); - printf ("(p in SI units): p_att = %g, p_MC= %g, p_ph1 = %g, p_ph2 = %g, p_ph3 = %g, p_J = %g, bose_factor = %g, V_rho = %g, J_factor/J0 = %g, G_factor = " - "%g +i %g, conj(G) = %g + I %g, |G|^2 = %g + I %g, ", - p_att, p_MC, p_ph1, p_ph2, p_ph3, p_J, bose_factor, V_rho, J_factor / J0, creal (Gn), cimag (Gn), creal (conj (Gn)), cimag (conj (Gn)), - creal (Gn * conj (Gn)), cimag (Gn * conj (Gn))); - printf (" p_tot = %g, p_final = %g\n", p_tot, p); - printf ("q is (h,k,l) = (%g, %g, %g) \n", q_x / astar, q_y / astar, q_z / cstar); - } - } /* else transmit: Neutron did not hit the sample */ + } %} MCDISPLAY From 93b54b8958852d45db8ae4be4041b07900733bbb Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Wed, 25 Feb 2026 22:10:17 +0100 Subject: [PATCH 15/15] Add %Example expectation value --- .../Tests_samples/Test_Phonon_BvK_PG/Test_Phonon_BvK_PG.instr | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/Test_Phonon_BvK_PG.instr b/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/Test_Phonon_BvK_PG.instr index aafff80bd..6cc880f46 100644 --- a/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/Test_Phonon_BvK_PG.instr +++ b/mcstas-comps/examples/Tests_samples/Test_Phonon_BvK_PG/Test_Phonon_BvK_PG.instr @@ -13,7 +13,7 @@ * %D * Simple TAS to test the Phonon_BvK_PG component * -* Example: +* %Example: -y Detector: e_mon_beforeana_zoom_I=6.90143e-16 * * %P * E: [meV] Energy transfer