Get mue

.gitmodules

@@ -1,3 +0,0 @@
-[submodule "NuLib"]
-	path = NuLib
-	url = git@github.com:srichers/NuLib.git

Makefile

@@ -0,0 +1,17 @@
+#!/bin/csh
+NULIB_LOC = /Users/jspark971/external/NuLib
+HDF5_LOC = /Users/jspark971/external/hdf5-1.10.5/hdf5
+# export DYLD_LIBRARY_PATH=/Users/jspark971/Documents/research/isotropicSQA/IsotropicSQA/external/hdf5-1.10.5/hdf5/lib/
+
+INCLUDE  = -I${HDF5_LOC}/include
+
+LIBRARY = ${NULIB_LOC}/src/nulib.a -L${HDF5_LOC}/lib -lhdf5 -lhdf5_fortran -lhdf5_cpp -L/usr/local/Cellar/gcc/9.2.0_3/lib/gcc/9 -lgfortran -lomp
+
+WHICHCODE = sqa2.psma
+
+COMP = g++ -g -std=gnu++11 -O3 -Wall -Xpreprocessor -fopenmp # -DNDEBUG
+
+${WHICHCODE}.o: ${WHICHCODE}.cpp
+	rm -f sqa2.psma.x sqa.tar.gz
+	tar --exclude=*.lum --exclude=*.cyl -cvzf sqa.tar.gz ./*
+	${COMP} ${WHICHCODE}.cpp -o ${WHICHCODE}.x ${INCLUDE} ${LIBRARY}

headers/interact.h

@@ -36,7 +36,7 @@ void initialize(array<array<MATRIX<complex<double>,NF,NF>,NE>,NM>& fmatrixf,
   cout << "T = " << T << " MeV" << endl;
   cout << "Ye = " << Ye << endl;
   eas.set(rho,T,Ye,do_interact);
-  
+
   for(int i=0; i<NE; i++){
     for(int m=matter; m<=antimatter; m++){
       double fe = eas.fermidirac(m, i);
@@ -47,7 +47,7 @@ void initialize(array<array<MATRIX<complex<double>,NF,NF>,NE>,NM>& fmatrixf,
       double z = (fe-fx)/2.;
       double xmax = sqrt(lmax*lmax - z*z);
       assert(xmax==xmax);
-      
+
       fmatrixf[m][i][e ][e ] = fe;
       fmatrixf[m][i][mu][mu] = fx;
       fmatrixf[m][i][e ][mu] = xmax*mixing;
@@ -147,7 +147,7 @@ array<array<MATRIX<complex<double>,NF,NF>,NE>,NM> my_interact
       for(flavour f1=e; f1<=mu; f1++)
 	for(flavour f2=e; f2<=mu; f2++)
 	  dfdr[m][i][f1][f2] = 0;
-      
+
       // emission
       dfdr[m][i][e ][e ] += eas.emis(se,i);
       dfdr[m][i][mu][mu] += eas.emis(sx,i);
@@ -157,16 +157,16 @@ array<array<MATRIX<complex<double>,NF,NF>,NE>,NM> my_interact
       for(flavour f1=e; f1<=mu; f1++)
 	for(flavour f2=e; f2<=mu; f2++)
 	  dfdr[m][i][f1][f2] -= kappa_abs[f1][f2] * fmatrixf[m][i][f1][f2];
-    
+
       // scattering and pair annihilation
       for(int j=0; j<NE; j++){
-	
+
 	// reusable variables
 	double Phi0e, Phi0x, conv_to_in_rate;
 	complex<double> unblock_in, unblock_out;
 	MATRIX<double,NF,NF> Phi0avg, Phi0tilde, Phi0;
 	MATRIX<complex<double>,NF,NF> block;
-      
+
 	// scattering rate from group i to group j
 	Phi0e = eas.Phi0scat(se,i,j); // cm^3/s/sr
 	Phi0x = eas.Phi0scat(sx,i,j);
@@ -184,7 +184,7 @@ array<array<MATRIX<complex<double>,NF,NF>,NE>,NM> my_interact
 	      (unblock_in - unblock_out - block[f1][f2]*(conv_to_in_rate - 1.));
 	  }
 	}
-      
+
 	// annihilation from group i and anti group j
 	if(eas.do_pair){
 	  Phi0e = eas.Phi0pair(se,i,j); // cm^3/s/sr
@@ -206,16 +206,16 @@ array<array<MATRIX<complex<double>,NF,NF>,NE>,NM> my_interact
 	    }
 	  }
 	}
-	
+
 	// 4-neutrino scattering
 	if(eas.do_nu4scat){
 	  for(int j3=0; j3<NE; j3++){
 	    int j2 = eas.nu4_bin2(i,j,j3);
 	    if(j2<0 or j2>=NE) continue;
-	    
-	    if(abs(fmatrixf[m][j][e][mu]) > 0)
-	      assert( abs(fmatrixf[m][j][e][mu] - conj(fmatrixf[m][j][mu][e])) / abs(fmatrixf[m][j][e][mu]) < 1e-5);
-	    
+
+	    if(abs(fmatrixf[m][j2][e][mu]) > 0) // j -> j2
+	      assert( abs(fmatrixf[m][j2][e][mu] - conj(fmatrixf[m][j2][mu][e])) / abs(fmatrixf[m][j2][e][mu]) < 1e-5);
+
 	    int index = eas.nu4_kernel_index(i,j,j3);
 	    double kernel = __nulibtable_MOD_nulibtable_nu4scat[index];
 	    index = eas.nu4_kernel_index(i,j3,j);
@@ -227,8 +227,15 @@ array<array<MATRIX<complex<double>,NF,NF>,NE>,NM> my_interact
 	    tmp = fmatrixf[m][j]*(1.-fmatrixf[m][j2]);
 	    Pi_plus  += (Trace(tmp) + tmp) * fmatrixf[m][j3] * kernel;
 	  }
+
+		if(abs(Pi_plus[e][mu]) > 0) {
+cout<<Pi_plus<< endl;
+assert( abs(fmatrixf[m][j][e][mu] - conj(fmatrixf[m][j][mu][e])) / abs(fmatrixf[m][j][e][mu]) < 1e-5);
+// assert( abs(fmatrixf[m][j2][e][mu] - conj(fmatrixf[m][j2][mu][e])) / abs(fmatrixf[m][j2][e][mu]) < 1e-5);
+assert( abs(Pi_plus[e][mu] - conj(Pi_plus[mu][e])) / abs(Pi_plus[e][mu]) < 1e-5);
+		}
 	}
-	
+
 	// 4-neutrino pair processes
 	if(eas.do_nu4pair){
 	  for(int j3=0; j3<NE; j3++){
@@ -253,13 +260,20 @@ array<array<MATRIX<complex<double>,NF,NF>,NE>,NM> my_interact
 	    tmp = fmatrixf[m][j3]*fmatrixf[mbar][j];
 	    Pi_plus += (Trace(tmp) + tmp) * (1.-fmatrixf[mbar][j2]) * kernel;
 	  }
+
+		if(abs(Pi_plus[e][mu]) > 0) {
+cout<<Pi_plus<< endl;
+assert( abs(Pi_plus[e][mu] - conj(Pi_plus[mu][e])) / abs(Pi_plus[e][mu]) < 1e-5);
+		}
 	}
-	  
+
 	assert(j<NE);
       }
 
-      if(abs(Pi_plus[e][mu]) > 0)
+      if(abs(Pi_plus[e][mu]) > 0) {
+	cout<<Pi_plus<< endl;
 	assert( abs(Pi_plus[e][mu] - conj(Pi_plus[mu][e])) / abs(Pi_plus[e][mu]) < 1e-5);
+			}
       if(abs(Pi_minus[e][mu]) > 0)
 	assert( abs(Pi_minus[e][mu] - conj(Pi_minus[mu][e])) / abs(Pi_minus[e][mu]) < 1e-5);
       dfdr[m][i] += Pi_plus *(1.-fmatrixf[m][i]) + (1.-fmatrixf[m][i])*Pi_plus ;

headers/nulib_interface.h

@@ -262,6 +262,19 @@ class EAS{
     free(__nulibtable_MOD_nulibtable_epannihiltable_phi1);
   }
 
+  double get_mue(const double rho /* g/ccm */, const double temp /*MeV*/, const double ye) const{ // MeV
+    double eos_variables[__nulib_MOD_total_eos_variables];
+    for(int i = 0; i<__nulib_MOD_total_eos_variables; i++)
+    eos_variables[i] = 0;
+    eos_variables[0] = rho;
+    eos_variables[1] = temp;
+    eos_variables[2] = ye;
+
+    set_eos_variables_(eos_variables);
+    double mue = eos_variables[10];
+    return mue;
+  }
+
   double get_munue(const double rho /* g/ccm */, const double temp /*MeV*/, const double ye) const{ // MeV
     double eos_variables[__nulib_MOD_total_eos_variables];
     for(int i=0; i<__nulib_MOD_total_eos_variables; i++) eos_variables[i] = 0;

makefiles/Makefile_Andy

@@ -1,6 +1,6 @@
 #!/bin/csh
-NULIB_LOC = external/NuLib
-HDF5_LOC = external/hdf5-1.10.5/hdf5
+NULIB_LOC = /Users/jspark971/external/NuLib
+HDF5_LOC = /Users/jspark971/external/hdf5-1.10.5/hdf5
 
 INCLUDE  = -I${HDF5_LOC}/include
 

plot/plotf.py

@@ -0,0 +1,135 @@
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+import sys
+import matplotlib as mpl
+import matplotlib.pylab as pylab
+from multiprocessing import pool
+from matplotlib.ticker import AutoMinorLocator
+import matplotlib.colors as colors
+import h5py
+
+def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=100):
+    new_cmap = colors.LinearSegmentedColormap.from_list(
+        'trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name, a=minval, b=maxval),
+        cmap(np.linspace(minval, maxval, n)))
+    return new_cmap
+orig_cmap = plt.get_cmap('gist_heat')
+cmap = truncate_colormap(orig_cmap, 0, 0.85, 200)
+
+clight = 3e10 # cm/s
+h = 6.6260755e-27 # erg s
+eV = 1.60218e-12 # erg
+
+ns = 2
+ng=25
+
+mpl.rcParams['font.size'] = 22
+mpl.rcParams['font.family'] = 'serif'
+mpl.rc('text', usetex=True)
+mpl.rcParams['xtick.major.size'] = 7
+mpl.rcParams['xtick.major.width'] = 2
+mpl.rcParams['xtick.major.pad'] = 8
+mpl.rcParams['xtick.minor.size'] = 4
+mpl.rcParams['xtick.minor.width'] = 2
+mpl.rcParams['ytick.major.size'] = 7
+mpl.rcParams['ytick.major.width'] = 2
+mpl.rcParams['ytick.minor.size'] = 4
+mpl.rcParams['ytick.minor.width'] = 2
+mpl.rcParams['axes.linewidth'] = 2
+
+E = (np.arange(0,25)+1)*4 #np.array([2.0, 2.43154, 2.95619, 3.59405, 4.36953, 5.31235, 6.45859, 7.85216, 9.54642, 11.6062, 14.1105, 17.1551, 20.8567, 25.3569, 30.8282, 37.48])
+
+def line_color(ig,ng, cmap):
+    return cmap(E[ig]/E[-1])
+
+def column(a,b, state, group, filename):
+    realindex = 1 + group*2*ns*ns*2 + state*ns*ns*2 + b*ns*2 + a*2
+    imagindex = realindex+1
+    realcol = np.genfromtxt(filename, usecols=(realindex))
+    imagcol = np.genfromtxt(filename, usecols=(imagindex))[:len(realcol)]
+    return realcol, imagcol
+
+def isospin(f):
+    fx = (f[:,:,0,1]+f[:,:,1,0])/2.
+    fz = (f[:,:,0,0]-f[:,:,1,1])/2.
+    ft = (f[:,:,0,0]+f[:,:,1,1])/2.
+    return fx,fz,ft
+
+def magnitude(r,i):
+    return np.sqrt(r**2 + i**2)
+
+def makeplot(location):
+    filename = location+"/f.h5"
+    f=h5py.File(filename,"r")
+    fig, axes = plt.subplots(4, 1, sharex=True, figsize=(8,12))
+    t = np.array(f["r(cm)"]) / clight * 1e12
+    fmatrixf = np.array(f["fmatrixf"])
+
+    for ig in range(ng):
+        r = fmatrixf[:,0,ig,0,1,0]
+        axes[0].plot(t, fmatrixf[:,0,ig,0,0,0], '-', color=line_color(ig,ng,cmap), linewidth=1.5)
+        p, = axes[1].plot(t, fmatrixf[:,0,ig,1,1,0], '-', color=line_color(ig,ng,cmap), linewidth=1.5)
+        axes[2].plot(t, fmatrixf[:,0,ig,1,0,0]/fmatrixf[0,0,ig,1,0,0], '-', color=line_color(ig,ng,cmap), linewidth=1.5)
+        axes[3].plot(t, fmatrixf[:,0,ig,1,0,1]/fmatrixf[0,0,ig,1,0,0], '-', color=line_color(ig,ng,cmap), linewidth=1.5)
+        if(ig==0):
+            p1=p
+    for ig in range(ng):
+        r = fmatrixf[:,1,ig,0,1,0]
+        axes[0].plot(t, fmatrixf[:,1,ig,0,0,0], '--', color=line_color(ig,ng,cmap), linewidth=1.5)
+        p, = axes[1].plot(t, fmatrixf[:,1,ig,1,1,0], '--', color=line_color(ig,ng,cmap), linewidth=1.5)
+        axes[2].plot(t, fmatrixf[:,1,ig,1,0,0]/fmatrixf[0,1,ig,1,0,0], '--', color=line_color(ig,ng,cmap), linewidth=1.5)
+        axes[3].plot(t, fmatrixf[:,1,ig,1,0,1]/fmatrixf[0,1,ig,1,0,0], '--', color=line_color(ig,ng,cmap), linewidth=1.5)
+        if ig==0:
+            p2=p
+            axes[1].legend([p1,p2],[r"$f$",r"$\bar{f}$"],fontsize=22, ncol=2,loc='upper center',bbox_to_anchor=(0.5, 1), handletextpad=0, frameon=False)
+
+    #axes[2].grid()
+    axes[3].set_xlabel(r"$t\,(\mathrm{ps})$",fontsize=22)
+    #axes[2].set_ylabel(r"$f_{e\mu}(t)/f_{e\mu}(0)$",fontsize=22)
+    #axes[2].savefig("foffdiag.pdf",bbox_inches="tight")
+
+
+    #####################################
+
+
+    #if ig==0:
+    #    axes[0].legend([p1,p2],[r"$f_{ee}$",r"$f_{\mu\mu}$"], fontsize=22,loc="upper right", ncol=2, frameon=False)
+    #    axes[1].legend([p3,p4],[r"$\bar{f}_{ee}$",r"$\bar{f}_{\mu\mu}$"],fontsize=22,loc="upper right", ncol=2, frameon=False)
+
+
+    # colorbar
+    cax = fig.add_axes([.92, 0.1, 0.04, 0.8])
+    norm = mpl.colors.Normalize(vmin=0, vmax=E[-1])
+    cb1 = mpl.colorbar.ColorbarBase(cax, cmap=cmap,norm=norm,orientation='vertical')
+    cb1.set_label(r'$h\nu\,(\mathrm{MeV})$')
+    cax.minorticks_on()
+
+    #axes[0].set_xlabel(r"$t\,\mathrm{(ms)}$",fontsize=22)
+    #axes[0].set_ylabel(r"$f(t)/f(0)$",fontsize=22)
+    ylabelx = -.15
+    plt.text(ylabelx,0.5,r"$f_{ee}(t)$",horizontalalignment='center',verticalalignment='center', transform=axes[0].transAxes,rotation=90,fontsize=22)
+    plt.text(ylabelx,0.5,r"$f_{\mu\mu}(t)$",horizontalalignment='center',verticalalignment='center', transform=axes[1].transAxes,rotation=90,fontsize=22)
+    plt.text(ylabelx,0.5,r"$f_{(x)}(t)/f_{(x)}(0)$",horizontalalignment='center',verticalalignment='center', transform=axes[2].transAxes,rotation=90,fontsize=22)
+    plt.text(ylabelx,0.5,r"$f_{(y)}(t)/f_{(x)}(0)$",horizontalalignment='center',verticalalignment='center', transform=axes[3].transAxes,rotation=90,fontsize=22)
+    plt.subplots_adjust(wspace=0, hspace=0)
+    plt.minorticks_on()
+    axes[0].yaxis.set_major_locator(plt.MultipleLocator(.5))
+    axes[1].yaxis.set_major_locator(plt.MultipleLocator(.2))
+    axes[2].yaxis.set_major_locator(plt.MultipleLocator(.5))
+    axes[3].yaxis.set_major_locator(plt.MultipleLocator(.5))
+    #axes[0].set_ylim(0,1.001)
+    #axes[1].set_ylim(0,1)
+    axes[2].set_ylim(-1.4,1.4)
+    axes[3].set_ylim(-1.4,1.4)
+    for ax in axes:
+        ax.xaxis.set_minor_locator(AutoMinorLocator())
+        ax.yaxis.set_minor_locator(AutoMinorLocator())
+        #ax.ticklabel_format(useOffset=False)
+        ax.set_xlim(0,.68)
+
+    plt.savefig(location+"/f_short.pdf",bbox_inches="tight")
+
+
+for loc in ["."]:#["noosc_abs","noosc_fakepair","noosc_pair","noosc_fakebrems","noosc_escat","noosc_iscat","noosc_nointeract","noosc_nucscat"]:
+    makeplot(loc)

sqa2.psma.cpp

@@ -109,13 +109,13 @@ int main(int argc, char *argv[]){
     output_file = setup_file(outputfilename, s);
     write_data(output_file, s, 0);
   }
-  
+
   // random number generator - prevent aliasing in interactions and output
   srand(time(NULL));
 
   // temporary variable
   array<array<MATRIX<complex<double>,NF,NF>,NE>,NM> fmatrixf0 = s.fmatrixf;
-  
+
 
   // ***********************
   // start the loop over r *
@@ -192,7 +192,7 @@ int main(int argc, char *argv[]){
     for(int m=matter;m<=antimatter;m++)
       for(int i=0;i<=s.eas.ng-1;i++)
 	Hermitize(s.fmatrixf[m][i], accuracy);
-    
+
   }while(finish==false);
 
   cout << endl << "Finished" << endl;