Bsweep trEPR with new defin of lw

from_bittl/src/p2.h

+#define	PI	M_PI
+#define DataT double
+#define MAX_NUMBER_POINTS 32
+#define BE 1.399624167     /* beta/h; (vorher 1.39968)  SZ 9/8/96 */
+#define GBETA (2.0023193*BE) 
+#define FACTOR 0.398942   /*(1/sqrt(2*PI))*/
+#define SQR(x)	((x)*(x))
+#define asc ".asc\0"
+#define deg (180.0/PI)
+#define rad (PI/180.0)
+/*#define register  short*/  /* SZ 26/8/95 */
+
+typedef DataT
+   vector[3],
+   vector4[4],
+   vector10[10];
+
+typedef vector
+   matrix[3];
+
+typedef char
+   strg[16];
+
+strg		filename[9], orient_name[3];
+
+double   	*spec, **Transpec, **tens, *orient_spec;
+
+double		*tensor;
+
+matrix      	aA1, aP, gP, gP2, gA1, gA12, gaP2, ga2, Dm,
+                rPQ, rtPQ, rCQ, rtCQ, rBC, rtBC, Ra,
+		Rz180, Rz120, Rz240;
+
+vector		rpA1, euler_PQ, euler_zP, euler_zQ, zP, zQ;
+                
+
+vector4		intens, trans;
+
+vector10	hyperfine_prob_P, hyperfine_prob_A1;
+
+DataT		fmw, h1, h2, dh, phi1, theta1, J, d, Jpd, Jmd, geffP,
+		geffA1, aeffP, aeffA1, hbP, hbA1, hbPA1, miA1, miP,
+		lorentz_hb, max_freq;
+
+int 		imiA1, imiP, nt, np, step_size, number_of_protons,
+		n_protons_p, NumOfPoints, delta_w_steps, num_freq_points,
+		alpha_step, astep, num_spec;
+
+char    	reference_syst, print_orient, print_Powder,
+		print_ESEEM_FT,  outfile[20];
+				
+    

from_bittl/src/p2.txt

+9700.                MIKROWELLEN-FREQUENZ
+3450. 3470. 1024      MAGNETFELD
+181 181                   GENAUIGKEIT
+2.00300 0.00000 0.00000
+0.00000 2.00260 0.00000
+0.00000 0.00000 2.00230 G-TENSOR P700+
+0.2                    LINIENBREITE P
+2.00620 0.00000 0.00000
+0.00000 2.00510 0.00000
+0.00000 0.00000 2.00220 G-TENSOR Qb-
+0.2                     LINIENBREITE Q
+Q                       REFERENZSYSTEM
+0.0  90.0    		PHI,THETA
+83.0 128.0 10.0   	ALPHA,BETA,GAMMA
++0.01 -1.7            J,D
+Y                       ORIENTIERUNGSSPEKTREN?
+Y                       PULVERSPEKTRUM?
+N                       ESEEM_FT?
+test2                   OUTPUT_FILENAME
+31.25 512                MAX.ESEEM-FT-FREQUENZ,FREQ-PUNKTE
+256                      HYPERFEIN_SCHRITTE
+0.                     LORENTZ_BREITE(MHz)
+3                       ANZAHL_DER_PROTONEN_AM_A1-
+-10.25 -1.472  0.000
+-1.472 -11.95  0.000
+0.000  0.000 -9.000
+0                       ANZAHL_DER_PROTONEN_AM_P+
+ 0.1  0.0  0.0
+ 0.0  0.1  0.0
+ 0.0  0.0  0.1
+
+10.25  1.472  0.000
+ 1.472 11.95  0.000
+ 0.000  0.000 9.000
+
+-4.5  0.0  0.0      
+ 0.0 -4.0  0.0      
+ 0.0  0.0 -3.5       
+
+ HFC Tensor A1- (9.4,12.8,9.0) 30 grd gedreht.

from_bittl/src/p2_utils.c

+void PrintMatrix(a,title)
+matrix	a;
+char	title[];
+{
+   int	i, j;
+
+   printf("\n         %s \n",title);
+   for (i=0;i<=2;i++) {
+      printf("     ");
+      for (j=0;j<=2;j++)
+	 printf("%9.5f%c",a[i][j],' ');
+      printf("\n");
+   }
+}
+
+void PrintVec(a,title)
+vector	a;
+char	title[];
+{
+   int	i;
+
+   printf("\n         %s \n",title);
+   for (i=0;i<=2;i++)
+      printf("%9.5f%c",a[i],' ');
+   printf("\n");
+}
+
+void PrintVec10(a,title)
+vector10  a;
+char	  title[];
+{
+   int	i;
+
+   printf("\n         %s \n",title);
+   for (i=0;i<=4;i++)
+       printf("%9.5f%c",a[i],' ');
+   printf("\n");
+   for (i=5;i<=9;i++)
+       printf("%9.5f%c",a[i],' ');
+   printf("\n");
+   printf("\n");
+}
+
+void CopyMatrix (a,b)
+matrix	a, b;
+{
+   short  i, j;
+
+   for (i=0; i<3; i++)
+      for (j=0; j<3; j++)
+	 a[i][j] = b[i][j];
+}
+
+void MulMat(a,b,c)
+matrix 	a, b, c;
+{
+   double  dummy;
+   matrix  d, e;
+   short i, j, k;
+
+   CopyMatrix(d,b);
+   CopyMatrix(e,c);
+   for (i=0; i<3; i++)
+      for (j=0; j<3; j++) {
+	 for (dummy=0, k=0; k<3; k++)
+	    dummy += d[i][k] * e[k][j];
+	 a[i][j] = dummy;
+      }
+}
+
+void Transpose (a,b)
+matrix	a, b;
+{
+   short i, j;
+   matrix  c;
+
+   CopyMatrix(c,b);
+   for (i=0; i<3; i++)
+      for (j=0; j<3; j++)
+	 a[j][i] = c[i][j];
+}
+
+void RotationMatrix(a,winkel,axis)
+matrix 	a;
+double	winkel;
+char	axis;
+{
+   short  i, j, k;
+
+   winkel = winkel * PI /180.0;
+   for (i=0; i<3; i++)
+      for (j=0; j<3; j++)
+	 a[i][j] = 0.0;
+   switch (axis) {
+   case 'x':
+      i=0;
+      j=1;
+      k=2;
+      break;
+   case 'y':
+      i=1;
+      j=2;
+      k=0;
+      break;
+   case 'z':
+      i=2;
+      j=0;
+      k=1;
+      break;
+   }
+   a[i][i] = 1.0;
+   a[j][j] = cos(winkel);
+   a[k][k] = cos(winkel);
+   a[j][k] = sin(winkel);
+   a[k][j] = -sin(winkel);
+}
+
+void NewLine(fp)
+FILE  *fp;
+{
+   char  c;
+
+   while ((c = getc(fp)) != '\n');
+}
+
+
+double Compliment(a)
+double	a;
+{
+   if (a*a > 1.0) {
+      printf(" Warning: Function \"Compliment\" called with a value greater than 1.0!! Returns a value of 0.0 ");
+      return 0.0;
+   }
+   else
+      return sqrt(1.0 - a*a);
+}
+
+void WriteSpec(low, high, sp, n,name)
+double    low, high, *sp;
+int	  n;
+char	  name[];
+{
+   int    i;
+   FILE   *fp, *fopen();
+   double range;
+
+   range = high - low;
+   strncat(name,asc,sizeof(asc));
+   if ((fp = fopen(name,"w")) == NULL)
+      printf("file %s cannot be opened",name);
+   for (i=0; i<n; i++)
+      fprintf(fp,"%.3f\t%.9f\n",low + i*range/(n-1), sp[i]);
+   fclose(fp);
+}
+

trEPR_spectrumIntegral.asv

+%
+clearvars
+
+Sys.S = [1/2 1/2];
+Sys.g = [2, 2.0075];
+
+Sys.initState = 'singlet';
+
+% Exchange coupling
+nj = 100;
+Js = linspace(1e-2, 100, nj);  % MHz
+
+gPlus = (Sys.g(1) + Sys.g(2))/2;
+mwFreq = 9.7;  % GHz
+B0 = planck/bmagn/gPlus*mwFreq*1e12;  % mT
+nPoint = 5000;
+msAx = linspace(9.4, 10, nPoint);  % GHz, mw sweep axis
+bsAx = mhz2mt(msAx*1e3, gPlus);  % mT, b field sweep
+
+% 
+% Frequency sweep
+%
+mhzLw = 10;
+Sys.lw = mhzLw;  % MHz
+ 
+Exp.mwRange = [min(msAx), max(msAx)];
+Exp.Field = B0;
+Exp.nPoints = nPoint;
+Exp.Harmonic = 0;
+
+msSim = zeros(nj, nPoint);
+
+for ij = 1:nj
+    Sys.J = Js(ij);
+    msSim(ij, :) = pepper(Sys, Exp);
+    % msSim(ij, :) = msSim(ij, :)/max(abs(msSim(ij, :)));
+end
+
+msTotIntegral = sum(msSim, 2)/nPoint;
+
+%
+% Field sweep
+%
+clear("Exp")
+
+Sys.lw = mhz2mt(mhzLw);  % mT
+
+Exp.Range = [min(bsAx), max(bsAx)];
+Exp.mwFreq = mwFreq;
+Exp.nPoints = nPoint;
+Exp.Harmonic = 0;
+
+bsSim = zeros(nj, nPoint);
+
+for ij = 1:nj
+    Sys.J = Js(ij);
+    bsSim(ij, :) = pepper(Sys, Exp);
+    % bsSim(ij, :) = bsSim(ij, :)/max(abs(bsSim(ij, :)));
+end
+
+bsTotIntegral = sum(bsSim, 2)/nPoint;
+
+clf
+plot(Js, msTotIntegral, 'o-', 'DisplayName', 'Mw sweep')
+yyaxis right
+plot(Js, bsTotIntegral, 'o-', 'DisplayName', 'B0 sweep')
+xlabel('J / MHz')
+yyaxis left
+ylabel('Integral of trEPR spectrum')
+legend()
+
+%%
+
+ij = 100;
+clf
+plot(msAx, msSim(ij, :))
+yyaxis right
+plot(msAx, flip(bsSim(ij, :)))
+
+%%
+clearvars
+
+Sys.S = [1/2 1/2];
+Sys.g = [2, 2.0075];
+
+Sys.initState = 'singlet';
+
+% Dipolar coupling
+ndip = 100;
+dips = linspace(1, 100, ndip);  % MHz
+
+gPlus = (Sys.g(1) + Sys.g(2))/2;
+mwFreq = 9.7;  % GHz
+B0 = planck/bmagn/gPlus*mwFreq*1e12;  % mT
+nPoint = 5000;
+msAx = linspace(9.4, 10, nPoint);  % GHz, mw sweep axis
+bsAx = mhz2mt(msAx*1e3, gPlus);  % mT, b field sweep
+
+% 
+% Frequency sweep
+%
+mhzLw = 10;
+Sys.lw = mhzLw;  % MHz
+ 
+Exp.mwRange = [min(msAx), max(msAx)];
+Exp.Field = B0;
+Exp.nPoints = nPoint;
+Exp.Harmonic = 0;
+
+msSim = zeros(ndip, nPoint);
+
+for ij = 1:ndip
+    Sys.dip = dips(ij);
+    msSim(ij, :) = pepper(Sys, Exp);
+    % msSim(ij, :) = msSim(ij, :)/max(abs(msSim(ij, :)));
+end
+
+msTotIntegral = sum(msSim, 2)/nPoint;
+
+%
+% Field sweep
+%
+clear("Exp")
+
+Sys.lw = mhz2mt(mhzLw);  % mT
+
+Exp.Range = [min(bsAx), max(bsAx)];
+Exp.mwFreq = mwFreq;
+Exp.nPoints = nPoint;
+Exp.Harmonic = 0;
+
+bsSim = zeros(ndip, nPoint);
+
+for ij = 1:ndip
+    Sys.dip = dips(ij);
+    bsSim(ij, :) = pepper(Sys, Exp);
+    % bsSim(ij, :) = bsSim(ij, :)/max(abs(bsSim(ij, :)));
+end
+
+bsTotIntegral = sum(bsSim, 2)/nPoint;
+
+clf
+plot(dips, msTotIntegral, 'o-', 'DisplayName', 'Mw sweep')
+yyaxis right
+plot(dips, bsTotIntegral, 'o-', 'DisplayName', 'B0 sweep')
+xlabel('J / MHz')
+yyaxis left
+ylabel('Integral of trEPR spectrum')
+legend()
+
+%%
+ij = 100;
+clf
+plot(msAx, msSim(ij, :))
+yyaxis right
+plot(msAx, flip(bsSim(ij, :)))
+
+%%
+clearvars
+
+Sys.S = [1/2 1/2];
+Sys.g = [2, 2.0075];
+Sys.dip = 5;  % MHz
+
+Sys.initState = 'singlet';
+
+gPlus = (Sys.g(1) + Sys.g(2))/2;
+mwFreq = 9.7;  % GHz
+B0 = planck/bmagn/gPlus*mwFreq*1e12;  % mT
+nPoint = 5000;
+msAx = linspace(9.4, 10, nPoint);  % GHz, mw sweep axis
+bsAx = mhz2mt(msAx*1e3, gPlus);  % mT, b field sweep
+
+% 
+% Frequency sweep
+%
+mhzLw = 10;
+Sys.lw = mhzLw;  % MHz
+ 
+Exp.mwRange = [min(msAx), max(msAx)];
+Exp.Field = B0;
+Exp.nPoints = nPoint;
+Exp.Harmonic = 0;
+
+Opt.GridSize = [91 0];  % 181 orientations
+Opt.GridSymmetry = 'C1';
+
+% Powder average
+msSim = pepper(Sys, Exp, Opt);
+
+clear("Opt")
+
+Exp.CrystalSymmetry = 1;
+Exp.MolFrame = [0 0 0]*pi/180;
+Exp.SampleFrame = [0 0 0]*pi/180; 
+
+Opt.separate = 'orientations';
+rho = deg2rad(linspace(0, 180, 181));              % 10 degree increments
+solidAngleWeight = sin(rho')/sum(sin(rho), 'all');
+
+% Rotation around x
+Exp.SampleRotation = {'x', rho};
+msSimRotX = pepper(Sys, Exp, Opt);
+msPowderRotX  = sum(solidAngleWeight.*msSimRotX);
+
+% Rotation around y
+Exp.SampleRotation = {'y', rho};
+msSimRotY = pepper(Sys, Exp, Opt);
+msPowderRotY  = sum(solidAngleWeight.*msSimRotY);
+
+%
+% Field sweep
+%
+clear("Exp")
+
+Sys.lw = mhz2mt(mhzLw);  % mT
+
+Exp.Range = [min(bsAx), max(bsAx)];
+Exp.mwFreq = mwFreq;
+Exp.nPoints = nPoint;
+Exp.Harmonic = 0;
+
+clear("Opt")
+Opt.GridSize = [91 0];  % 181 orientations
+Opt.GridSymmetry = 'C1';
+
+% Powder average
+bsSim = pepper(Sys, Exp, Opt);
+
+clear("Opt")
+
+Exp.CrystalSymmetry = 1;
+Exp.MolFrame = [0 0 0]*pi/180;
+Exp.SampleFrame = [0 0 0]*pi/180; 
+
+Opt.separate = 'orientations';
+rho = deg2rad(linspace(0, 180, 181));
+solidAngleWeight = sin(rho')/sum(sin(rho), 'all');
+
+% Rotation around x
+Exp.SampleRotation = {'x', rho};
+bsSimRotX = pepper(Sys, Exp, Opt);
+bsPowderRotX  = sum(solidAngleWeight.*bsSimRotX);
+
+% Rotation around y
+Exp.SampleRotation = {'y', rho};
+bsSimRotY = pepper(Sys, Exp, Opt);
+bsPowderRotY  = sum(solidAngleWeight.*bsSimRotY);
+
+%%
+
+clf
+plot(msAx, msPowderRotX, msAx, msPowderRotY)
+% clf
+% plot(bsAx, bsPowderRotX, bsAx, bsPowderRotY)
+

trEPR_spectrumIntegral.m

+%% 1. Exhange coupling
+clearvars
+
+Sys.S = [1/2 1/2];
+Sys.g = [2, 2.0075];
+
+Sys.initState = 'singlet';
+
+% Exchange coupling
+nj = 100;
+Js = linspace(1e-2, 100, nj);  % MHz
+
+gPlus = (Sys.g(1) + Sys.g(2))/2;
+mwFreq = 9.7;  % GHz
+B0 = planck/bmagn/gPlus*mwFreq*1e12;  % mT
+nPoint = 5000;
+msAx = linspace(9.4, 10, nPoint);  % GHz, mw sweep axis
+bsAx = mhz2mt(msAx*1e3, gPlus);  % mT, b field sweep
+
+% 
+% Frequency sweep
+%
+mhzLw = 10;
+Sys.lw = mhzLw;  % MHz
+ 
+Exp.mwRange = [min(msAx), max(msAx)];
+Exp.Field = B0;
+Exp.nPoints = nPoint;
+Exp.Harmonic = 0;
+
+msSim = zeros(nj, nPoint);
+
+for ij = 1:nj
+    Sys.J = Js(ij);
+    msSim(ij, :) = pepper(Sys, Exp);
+    % msSim(ij, :) = msSim(ij, :)/max(abs(msSim(ij, :)));
+end
+
+msTotIntegral = sum(msSim, 2)/nPoint;
+
+%
+% Field sweep
+%
+clear("Exp")
+
+Sys.lw = mhz2mt(mhzLw);  % mT
+
+Exp.Range = [min(bsAx), max(bsAx)];
+Exp.mwFreq = mwFreq;
+Exp.nPoints = nPoint;
+Exp.Harmonic = 0;
+
+bsSim = zeros(nj, nPoint);
+
+for ij = 1:nj
+    Sys.J = Js(ij);
+    bsSim(ij, :) = pepper(Sys, Exp);
+    % bsSim(ij, :) = bsSim(ij, :)/max(abs(bsSim(ij, :)));
+end
+
+bsTotIntegral = sum(bsSim, 2)/nPoint;
+
+clf
+plot(Js, msTotIntegral, 'o-', 'DisplayName', 'Mw sweep')
+yyaxis right
+plot(Js, bsTotIntegral, 'o-', 'DisplayName', 'B0 sweep')
+xlabel('J / MHz')
+yyaxis left
+ylabel('Integral of trEPR spectrum')
+legend()
+
+%% Plot
+
+ij = 100;
+clf
+plot(msAx, msSim(ij, :))
+yyaxis right
+plot(msAx, flip(bsSim(ij, :)))
+
+%% Dipolar coupling
+clearvars
+
+Sys.S = [1/2 1/2];
+Sys.g = [2, 2.0075];
+
+Sys.initState = 'singlet';
+
+% Dipolar coupling
+ndip = 100;
+dips = linspace(1, 100, ndip);  % MHz
+
+gPlus = (Sys.g(1) + Sys.g(2))/2;
+mwFreq = 9.7;  % GHz
+B0 = planck/bmagn/gPlus*mwFreq*1e12;  % mT
+nPoint = 5000;
+msAx = linspace(9.4, 10, nPoint);  % GHz, mw sweep axis
+bsAx = mhz2mt(msAx*1e3, gPlus);  % mT, b field sweep
+
+% 
+% Frequency sweep
+%
+mhzLw = 10;
+Sys.lw = mhzLw;  % MHz
+ 
+Exp.mwRange = [min(msAx), max(msAx)];
+Exp.Field = B0;
+Exp.nPoints = nPoint;
+Exp.Harmonic = 0;
+
+msSim = zeros(ndip, nPoint);
+
+for ij = 1:ndip
+    Sys.dip = dips(ij);
+    msSim(ij, :) = pepper(Sys, Exp);
+    % msSim(ij, :) = msSim(ij, :)/max(abs(msSim(ij, :)));
+end
+
+msTotIntegral = sum(msSim, 2)/nPoint;
+
+%
+% Field sweep
+%
+clear("Exp")
+
+Sys.lw = mhz2mt(mhzLw);  % mT
+
+Exp.Range = [min(bsAx), max(bsAx)];
+Exp.mwFreq = mwFreq;
+Exp.nPoints = nPoint;
+Exp.Harmonic = 0;
+
+bsSim = zeros(ndip, nPoint);
+
+for ij = 1:ndip
+    Sys.dip = dips(ij);
+    bsSim(ij, :) = pepper(Sys, Exp);
+    % bsSim(ij, :) = bsSim(ij, :)/max(abs(bsSim(ij, :)));
+end
+
+bsTotIntegral = sum(bsSim, 2)/nPoint;
+
+clf
+plot(dips, msTotIntegral, 'o-', 'DisplayName', 'Mw sweep')
+yyaxis right
+plot(dips, bsTotIntegral, 'o-', 'DisplayName', 'B0 sweep')
+xlabel('J / MHz')
+yyaxis left
+ylabel('Integral of trEPR spectrum')
+legend()
+
+%% Plot
+
+ij = 100;
+clf
+plot(msAx, msSim(ij, :))
+yyaxis right
+plot(msAx, flip(bsSim(ij, :)))
+
+%% Crystal rotation
+clearvars
+
+Sys.S = [1/2 1/2];
+Sys.g = [2, 2.0075];
+Sys.dip = 5;  % MHz
+
+Sys.initState = 'singlet';
+
+gPlus = (Sys.g(1) + Sys.g(2))/2;
+mwFreq = 9.7;  % GHz
+B0 = planck/bmagn/gPlus*mwFreq*1e12;  % mT
+nPoint = 5000;
+msAx = linspace(9.4, 10, nPoint);  % GHz, mw sweep axis
+bsAx = mhz2mt(msAx*1e3, gPlus);  % mT, b field sweep
+
+% 
+% Frequency sweep
+%
+mhzLw = 10;
+Sys.lw = mhzLw;  % MHz
+ 
+Exp.mwRange = [min(msAx), max(msAx)];
+Exp.Field = B0;
+Exp.nPoints = nPoint;
+Exp.Harmonic = 0;
+
+Opt.GridSize = [91 0];  % 181 orientations
+Opt.GridSymmetry = 'C1';
+
+% Powder average
+msSim = pepper(Sys, Exp, Opt);
+
+clear("Opt")
+
+Exp.CrystalSymmetry = 1;
+Exp.MolFrame = [0 0 0]*pi/180;
+Exp.SampleFrame = [0 0 0]*pi/180; 
+
+Opt.separate = 'orientations';
+rho = deg2rad(linspace(0, 180, 181));              % 10 degree increments
+solidAngleWeight = sin(rho')/sum(sin(rho), 'all');
+
+% Rotation around x
+Exp.SampleRotation = {'x', rho};
+msSimRotX = pepper(Sys, Exp, Opt);
+msPowderRotX  = sum(solidAngleWeight.*msSimRotX);
+
+% Rotation around y
+Exp.SampleRotation = {'y', rho};
+msSimRotY = pepper(Sys, Exp, Opt);
+msPowderRotY  = sum(solidAngleWeight.*msSimRotY);
+
+%
+% Field sweep
+%
+clear("Exp")
+
+Sys.lw = mhz2mt(mhzLw);  % mT
+
+Exp.Range = [min(bsAx), max(bsAx)];
+Exp.mwFreq = mwFreq;
+Exp.nPoints = nPoint;
+Exp.Harmonic = 0;
+
+clear("Opt")
+Opt.GridSize = [91 0];  % 181 orientations
+Opt.GridSymmetry = 'C1';
+
+% Powder average
+bsSim = pepper(Sys, Exp, Opt);
+
+clear("Opt")
+
+Exp.CrystalSymmetry = 1;
+Exp.MolFrame = [0 0 0]*pi/180;
+Exp.SampleFrame = [0 0 0]*pi/180; 
+
+Opt.separate = 'orientations';
+rho = deg2rad(linspace(0, 180, 181));
+solidAngleWeight = sin(rho')/sum(sin(rho), 'all');
+
+% Rotation around x
+Exp.SampleRotation = {'x', rho};
+bsSimRotX = pepper(Sys, Exp, Opt);
+bsPowderRotX  = sum(solidAngleWeight.*bsSimRotX);
+
+% Rotation around y
+Exp.SampleRotation = {'y', rho};
+bsSimRotY = pepper(Sys, Exp, Opt);
+bsPowderRotY  = sum(solidAngleWeight.*bsSimRotY);
+
+%% Plot
+
+clf
+plot(msAx, msPowderRotX, msAx, msPowderRotY)
+% clf
+% plot(bsAx, bsPowderRotX, bsAx, bsPowderRotY)
+