diff --git a/gecatsim/clib_build/src/analytic_projector.c b/gecatsim/clib_build/src/analytic_projector.c
index bb9b5ad..4ccd5d4 100644
--- a/gecatsim/clib_build/src/analytic_projector.c
+++ b/gecatsim/clib_build/src/analytic_projector.c
@@ -32,7 +32,7 @@ int n_col_oversample_add_xtalk = 1;
 int n_row_oversample_add_xtalk = 1;
 
 
-struct module_info
+struct analytic_module_info
 {
   double *Height;
   double *Width;
@@ -48,7 +48,7 @@ struct module_info
   int moduleOverlapType;
 };
 
-struct module_info modules = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0,0,0,0};
+struct analytic_module_info analytic_modules = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0,0,0,0};
 
 struct pnt
 {
@@ -67,7 +67,7 @@ struct bounding_info
 
 struct bounding_info bounding = {NULL,NULL};
 
-struct phantom_info
+struct analytic_phantom_info
 {
   int numObjects;
   int *objectType;
@@ -83,16 +83,16 @@ struct phantom_info
   double xbounds[2];
 };
 
-struct phantom_info phantom = {0,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,{0,0}};
+struct analytic_phantom_info analytic_phantom = {0,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,{0,0}};
 
-struct material_info
+struct analytic_material_info
 {
   int materialCount;
   int eBinCount;
   double *muTable;
 };
 
-struct material_info materials = {0,0,NULL};
+struct analytic_material_info analytic_materials = {0,0,NULL};
 
 struct projector_args
 {
@@ -463,12 +463,12 @@ int quartic_intersect_C(double *a0,double *Qlp,double *Qrp,double shp,double *al
   int i,out;
   double scale,displ,a0t[3],tmp[3],aa,b,c,d,e,f,C[5],*Ql,*Qr,sh;
   
-  Ql = &phantom.Qmatrix[obj*18];
-  //  printf("%1.12lf %1.12lf",phantom.Qmatrix[obj);
+  Ql = &analytic_phantom.Qmatrix[obj*18];
+  //  printf("%1.12lf %1.12lf",analytic_phantom.Qmatrix[obj);
   //  for(i = 1;i<9;i++) printf("L %1.12lf %1.12lf diff: %1.12lf\r\n",Ql[i],Qlp[i],Ql[i]-Qlp[i]);
-  Qr = &phantom.Qmatrix[obj*18+9];
+  Qr = &analytic_phantom.Qmatrix[obj*18+9];
   //  for(i = 1;i<9;i++) printf("R %1.12lf %1.12lf diff: %1.12lf\r\n",Qr[i],Qrp[i],Qr[i]-Qrp[i]);
-  sh = phantom.shape[obj];
+  sh = analytic_phantom.shape[obj];
   //  printf("s diff: %1.12lf\r\n",sh-shp);
   
   scale = 1.0;
@@ -505,9 +505,9 @@ int quartic_intersect(double *a0,double *alpha,double *tc,int obj)
   int i,out;
   double scale,displ,a0t[3],tmp[3],aa,b,c,d,e,f,C[5],*Ql,*Qr,sh;
   
-  Ql = &phantom.Qmatrix[obj*18];
-  Qr = &phantom.Qmatrix[obj*18+9];
-  sh = phantom.shape[obj];
+  Ql = &analytic_phantom.Qmatrix[obj*18];
+  Qr = &analytic_phantom.Qmatrix[obj*18+9];
+  sh = analytic_phantom.shape[obj];
   scale = 1.0;
   displ = 0;
   for(i = 0;i<3;i++) displ += a0[i]*a0[i];
@@ -541,12 +541,12 @@ DLLEXPORT int clip_all(double *a,double *alpha,double rayLength,double *tc2,int
   double b[3], s1, s2, tcrit, tmin, tmax, *eta, *s, den;
   int j, cp, firstPlane, materialIndex;
 
-  firstPlane = phantom.clipStartIndex[i];
-  eta = &phantom.clipNormalVector[firstPlane*3];
-  s = &phantom.clipDistance[firstPlane];
-  cp = phantom.nClipPlanes[i];
-  den = phantom.density[i];
-  materialIndex = phantom.materialInd[i];
+  firstPlane = analytic_phantom.clipStartIndex[i];
+  eta = &analytic_phantom.clipNormalVector[firstPlane*3];
+  s = &analytic_phantom.clipDistance[firstPlane];
+  cp = analytic_phantom.nClipPlanes[i];
+  den = analytic_phantom.density[i];
+  materialIndex = analytic_phantom.materialInd[i];
   tmin = -VERY_BIG;tmax = VERY_BIG;
   for(j = 0;j<3;j++) b[j] = a[j]+rayLength*alpha[j];
   
@@ -604,8 +604,8 @@ DLLEXPORT int quadratic_intersect(double *a0,double *alpha,int pars11,double *tc
   int out;
   double A,B,C,tmp,*Q,k;
     // The quadratic formula is of the form    A t^2 + B t + C = 0
-  Q = &phantom.Qmatrix[obj*18];
-  k = phantom.shape[obj];
+  Q = &analytic_phantom.Qmatrix[obj*18];
+  k = analytic_phantom.shape[obj];
   C = quadratic_form(a0,Q,a0)-k;
   B = 2.0*quadratic_form(alpha,Q,a0);
   A = quadratic_form(alpha,Q,alpha);
@@ -665,29 +665,29 @@ double* my_memcpyd(double *src, double *dest, int bytes)
 DLLEXPORT void set_src_info(double *sourceWeights, int nSubSources)
 
 {
-  modules.sourceWeights = my_memcpyd(sourceWeights, modules.sourceWeights, nSubSources*sizeof(double));
-  modules.nSubSources = nSubSources;
+  analytic_modules.sourceWeights = my_memcpyd(sourceWeights, analytic_modules.sourceWeights, nSubSources*sizeof(double));
+  analytic_modules.nSubSources = nSubSources;
 }
 
 DLLEXPORT void set_module_info(double *Height, double *Width, int *Pix, double *Coords, int *Sub, double *Sampling, double *Weight, int nModuleTypes, int maxPix, int maxSubDets, int moduleOverlapType)
 
 {
   int i;
-  modules.Height = my_memcpyd(Height, modules.Height, nModuleTypes*sizeof(double));
-  modules.Width = my_memcpyd(Width, modules.Width, nModuleTypes*sizeof(double));
+  analytic_modules.Height = my_memcpyd(Height, analytic_modules.Height, nModuleTypes*sizeof(double));
+  analytic_modules.Width = my_memcpyd(Width, analytic_modules.Width, nModuleTypes*sizeof(double));
   for(i = 0; i < nModuleTypes; i++)
     {
-      modules.Height[i] = MAX(1e-7,modules.Height[i]);
-      modules.Width[i] = MAX(1e-7,modules.Width[i]);
+      analytic_modules.Height[i] = MAX(1e-7,analytic_modules.Height[i]);
+      analytic_modules.Width[i] = MAX(1e-7,analytic_modules.Width[i]);
     }
-  modules.Pix = my_memcpyi(Pix, modules.Pix, nModuleTypes*sizeof(int));
-  modules.Coords = my_memcpyd(Coords, modules.Coords, maxPix*2*nModuleTypes*sizeof(double));
-  modules.Sub = my_memcpyi(Sub, modules.Sub, nModuleTypes*sizeof(int));
-  modules.Sampling = my_memcpyd(Sampling, modules.Sampling, 2*maxSubDets*nModuleTypes*sizeof(double));
-  modules.Weight = my_memcpyd(Weight, modules.Weight, maxSubDets*nModuleTypes*sizeof(double));
-  modules.maxPixPerModule = maxPix;
-  modules.maxSubDets = maxSubDets;
-  modules.moduleOverlapType = moduleOverlapType;
+  analytic_modules.Pix = my_memcpyi(Pix, analytic_modules.Pix, nModuleTypes*sizeof(int));
+  analytic_modules.Coords = my_memcpyd(Coords, analytic_modules.Coords, maxPix*2*nModuleTypes*sizeof(double));
+  analytic_modules.Sub = my_memcpyi(Sub, analytic_modules.Sub, nModuleTypes*sizeof(int));
+  analytic_modules.Sampling = my_memcpyd(Sampling, analytic_modules.Sampling, 2*maxSubDets*nModuleTypes*sizeof(double));
+  analytic_modules.Weight = my_memcpyd(Weight, analytic_modules.Weight, maxSubDets*nModuleTypes*sizeof(double));
+  analytic_modules.maxPixPerModule = maxPix;
+  analytic_modules.maxSubDets = maxSubDets;
+  analytic_modules.moduleOverlapType = moduleOverlapType;
 }
 
 DLLEXPORT void set_bounding_info(int numObjs, int *vertexStartInd, double *vertLocs, int numVerts)
@@ -697,38 +697,38 @@ DLLEXPORT void set_bounding_info(int numObjs, int *vertexStartInd, double *vertL
 
   bounding.vertexStartIndex = my_memcpyi(vertexStartInd, bounding.vertexStartIndex, (numObjs+1)*sizeof(int));
   bounding.vertexLocations = my_memcpyd(vertLocs, bounding.vertexLocations, sizeof(double)*numVerts*3);
-  phantom.xbounds[0]=VERY_BIG;
-  phantom.xbounds[1]=-VERY_BIG;
+  analytic_phantom.xbounds[0]=VERY_BIG;
+  analytic_phantom.xbounds[1]=-VERY_BIG;
   for(i=0;i<numVerts;i++) {
-    if (bounding.vertexLocations[i*3]>phantom.xbounds[1]) phantom.xbounds[1]=bounding.vertexLocations[i*3];
-    if (bounding.vertexLocations[i*3]<phantom.xbounds[0]) phantom.xbounds[0]=bounding.vertexLocations[i*3];
+    if (bounding.vertexLocations[i*3]>analytic_phantom.xbounds[1]) analytic_phantom.xbounds[1]=bounding.vertexLocations[i*3];
+    if (bounding.vertexLocations[i*3]<analytic_phantom.xbounds[0]) analytic_phantom.xbounds[0]=bounding.vertexLocations[i*3];
   }
 }
 
 DLLEXPORT void set_phantom_info(int numObjs, int *objType, int *clipStInd, int *nPlanes, int *matInd, double *objCent, double *shp, double *Qmat, double *clipNormVec, double *clipDist, double *den, int totalNumPlanes)
 
 {
-  phantom.numObjects = numObjs;
-
-  phantom.objectType = my_memcpyi(objType, phantom.objectType, numObjs*sizeof(int));
-  phantom.clipStartIndex = my_memcpyi(clipStInd, phantom.clipStartIndex, (numObjs+1)*sizeof(int));
-  phantom.nClipPlanes = my_memcpyi(nPlanes, phantom.nClipPlanes, numObjs*sizeof(int));
-  phantom.materialInd = my_memcpyi(matInd, phantom.materialInd, numObjs*sizeof(int));
-
-  phantom.objectCenter = my_memcpyd(objCent, phantom.objectCenter, numObjs*sizeof(double)*3);
-  phantom.shape = my_memcpyd(shp, phantom.shape, numObjs*sizeof(double));
-  phantom.Qmatrix = my_memcpyd(Qmat, phantom.Qmatrix, numObjs*sizeof(double)*18);
-  phantom.clipNormalVector = my_memcpyd(clipNormVec, phantom.clipNormalVector, sizeof(double)*totalNumPlanes*3);
-  phantom.clipDistance = my_memcpyd(clipDist, phantom.clipDistance, sizeof(double)*totalNumPlanes);
-  phantom.density = my_memcpyd(den, phantom.density, numObjs*sizeof(double));
+  analytic_phantom.numObjects = numObjs;
+
+  analytic_phantom.objectType = my_memcpyi(objType, analytic_phantom.objectType, numObjs*sizeof(int));
+  analytic_phantom.clipStartIndex = my_memcpyi(clipStInd, analytic_phantom.clipStartIndex, (numObjs+1)*sizeof(int));
+  analytic_phantom.nClipPlanes = my_memcpyi(nPlanes, analytic_phantom.nClipPlanes, numObjs*sizeof(int));
+  analytic_phantom.materialInd = my_memcpyi(matInd, analytic_phantom.materialInd, numObjs*sizeof(int));
+
+  analytic_phantom.objectCenter = my_memcpyd(objCent, analytic_phantom.objectCenter, numObjs*sizeof(double)*3);
+  analytic_phantom.shape = my_memcpyd(shp, analytic_phantom.shape, numObjs*sizeof(double));
+  analytic_phantom.Qmatrix = my_memcpyd(Qmat, analytic_phantom.Qmatrix, numObjs*sizeof(double)*18);
+  analytic_phantom.clipNormalVector = my_memcpyd(clipNormVec, analytic_phantom.clipNormalVector, sizeof(double)*totalNumPlanes*3);
+  analytic_phantom.clipDistance = my_memcpyd(clipDist, analytic_phantom.clipDistance, sizeof(double)*totalNumPlanes);
+  analytic_phantom.density = my_memcpyd(den, analytic_phantom.density, numObjs*sizeof(double));
 }
 
 DLLEXPORT void set_material_info(int materialCount, int eBinCount, double *muTable)
 
 {
-  materials.materialCount = materialCount;
-  materials.eBinCount = eBinCount;
-  materials.muTable = my_memcpyd(muTable, materials.muTable, eBinCount*materialCount*sizeof(double));
+  analytic_materials.materialCount = materialCount;
+  analytic_materials.eBinCount = eBinCount;
+  analytic_materials.muTable = my_memcpyd(muTable, analytic_materials.muTable, eBinCount*materialCount*sizeof(double));
 }
 
 DLLEXPORT void intersections_full_list(int *objlist, int lenObjList,double *a,double *alpha,double rayLength,double *t_ends,int *matls,double *dens,int *num_segs)
@@ -743,8 +743,8 @@ DLLEXPORT void intersections_full_list(int *objlist, int lenObjList,double *a,do
 
   for(n = 0;n<lenObjList;n++){
     i = objlist[n];         // Loop over objects of interest
-    p11 = phantom.objectType[i];
-    for(j = 0;j<3;j++) a0[j] = (a[j]-phantom.objectCenter[i*3+j]);              // vector pointing from obj. center to src
+    p11 = analytic_phantom.objectType[i];
+    for(j = 0;j<3;j++) a0[j] = (a[j]-analytic_phantom.objectCenter[i*3+j]);              // vector pointing from obj. center to src
     for(j = 0;j<4;j++) tc2[j] = 0.0;
     if((p11 != 3) && (p11 != 7))                  // If not a torus or vessel segment
       out = quadratic_intersect(a0,alpha,p11,tc2,i);
@@ -875,30 +875,30 @@ DLLEXPORT void make_vol(float *volume, int Nx, double xoff, double dx, int Ny, d
   int i, j, k, l, m, p, *matls, *objlist, num_segs, energyBinsToSkip = 0, volume_offset;
   double a[3], alpha[3], *t_ends, *dens, rayLength, xbrd, *matl_tabl;
 
-  matl_tabl = malloc(sizeof(double)*materials.materialCount);
-  matls = malloc(sizeof(int)*phantom.numObjects*4);
-  objlist = malloc(sizeof(int)*phantom.numObjects);
-  t_ends = malloc(sizeof(double)*phantom.numObjects*4);
-  dens = malloc(sizeof(double)*phantom.numObjects*4);
-  for(i = 0;i<materials.materialCount;i++) {
+  matl_tabl = malloc(sizeof(double)*analytic_materials.materialCount);
+  matls = malloc(sizeof(int)*analytic_phantom.numObjects*4);
+  objlist = malloc(sizeof(int)*analytic_phantom.numObjects);
+  t_ends = malloc(sizeof(double)*analytic_phantom.numObjects*4);
+  dens = malloc(sizeof(double)*analytic_phantom.numObjects*4);
+  for(i = 0;i<analytic_materials.materialCount;i++) {
     if (material_volumes) 
       matl_tabl[i] = 1.0/(oversampling*oversampling);
     else
-      matl_tabl[i] = materials.muTable[i+materials.materialCount*energyBinsToSkip]/(oversampling*oversampling);
+      matl_tabl[i] = analytic_materials.muTable[i+analytic_materials.materialCount*energyBinsToSkip]/(oversampling*oversampling);
     printf("\n\rmat%d:  %1.12lf\n\r",i,matl_tabl[i]*oversampling*oversampling);
   }
   dy = dy/oversampling;
   yoff = oversampling*(yoff-1)+(oversampling+1)*0.5;
   dz = dz/oversampling;
   zoff = oversampling*(zoff-1)+(oversampling+1)*0.5;
-  for(i = 0;i<phantom.numObjects;i++) objlist[i] = i;
+  for(i = 0;i<analytic_phantom.numObjects;i++) objlist[i] = i;
 
   alpha[0] = 1;
   alpha[1] = 0;
   alpha[2] = 0;
-  a[0] = phantom.xbounds[0] - 1;
-  rayLength = phantom.xbounds[1] + 1 - phantom.xbounds[0];
-  dbug(1,"\n\n\rphantom.xbounds[0]:  %1.12lf  \n\r phantom.xbounds[0]:  %1.12lf  \n\rmaterial_volumes: %d",phantom.xbounds[0],phantom.xbounds[1],material_volumes);
+  a[0] = analytic_phantom.xbounds[0] - 1;
+  rayLength = analytic_phantom.xbounds[1] + 1 - analytic_phantom.xbounds[0];
+  dbug(1,"\n\n\ranalytic_phantom.xbounds[0]:  %1.12lf  \n\r analytic_phantom.xbounds[0]:  %1.12lf  \n\rmaterial_volumes: %d",analytic_phantom.xbounds[0],analytic_phantom.xbounds[1],material_volumes);
   for(l = 0;l<Nz;l++)
     for(k = l*oversampling;k<l*oversampling+oversampling;k++){
       a[2] = (k+1-zoff)*dz;
@@ -906,7 +906,7 @@ DLLEXPORT void make_vol(float *volume, int Nx, double xoff, double dx, int Ny, d
 	for(j = m*oversampling;j<m*oversampling+oversampling;j++){
 	  a[1] = (j+1-yoff)*dy;
 	  num_segs = 0;
-	  intersections_full_list(objlist, phantom.numObjects, a, alpha, rayLength, t_ends, matls, dens, &num_segs);
+	  intersections_full_list(objlist, analytic_phantom.numObjects, a, alpha, rayLength, t_ends, matls, dens, &num_segs);
 	  for(i = 0;i<num_segs;i++){
 	    t_ends[i] = t_ends[i]+a[0];
 	    //printf("%1.12lf\r\n",t_ends[i]);
@@ -961,8 +961,8 @@ DLLEXPORT void intersections(int *objlist,int lenObjList,double *a,double *alpha
   for(n = 0;n<lenObjList;n++){
     i = objlist[n];         // Loop over objects of interest
     //    if((i == 4) && (prnt)) pnt = 1; else pnt = 0;
-    p11 = phantom.objectType[i];
-    for(j = 0;j<3;j++) a0[j] = (a[j]-phantom.objectCenter[i*3+j]);              // vector pointing from obj. center to src
+    p11 = analytic_phantom.objectType[i];
+    for(j = 0;j<3;j++) a0[j] = (a[j]-analytic_phantom.objectCenter[i*3+j]);              // vector pointing from obj. center to src
     for(j = 0;j<4;j++) tc2[j] = 0.0;
     if((p11 != 3) && (p11 != 7))                  // If not a torus or vessel segment
       out = quadratic_intersect(a0,alpha,p11,tc2,i);
@@ -1644,7 +1644,7 @@ void store_half_planes(double *pixel_R, double *pixel_U, int *indexList, int lis
 void store(double *pixel_R, double *pixel_U, int *indexList, int listLength, int objectNumber, void *boundaries)
 
 {
-  switch(modules.moduleOverlapType)
+  switch(analytic_modules.moduleOverlapType)
     {
     case 1:
       store_height_lims(pixel_U, indexList, listLength, objectNumber, boundaries);
@@ -1679,12 +1679,12 @@ void compute_object_projections(double *srcHullPoints, int nSrcHullPoints, int m
   int *vertexPoints = NULL;
   int planB;
 
-  vertexPoints = malloc(sizeof(int)*phantom.numObjects);
+  vertexPoints = malloc(sizeof(int)*analytic_phantom.numObjects);
   //dbug(3,"COP1\r\n");
   cross(right,up,normal);
 
   // Compute the maximum number of vertex points on a bounding polyhedron
-  for(i = 0;i<phantom.numObjects;i++)
+  for(i = 0;i<analytic_phantom.numObjects;i++)
     {
       vertexPoints[i] = bounding.vertexStartIndex[i+1] - bounding.vertexStartIndex[i];
       if (vertexPoints[i]>maxVertexPoints) maxVertexPoints = vertexPoints[i];
@@ -1695,7 +1695,7 @@ void compute_object_projections(double *srcHullPoints, int nSrcHullPoints, int m
   indexList = malloc(sizeof(int)*(maxVertexPoints*nSrcHullPoints+moduleEdges*2));
 
   //dbug(3,"COP1\r\n");
-  for(i = 0;i<phantom.numObjects;i++)
+  for(i = 0;i<analytic_phantom.numObjects;i++)
     {
       //dbug(3,"Loop %d\r\n",i);
       planB = 0;
@@ -1763,7 +1763,7 @@ void compute_object_projections(double *srcHullPoints, int nSrcHullPoints, int m
       //dbug(3,"loopend2 %d %d\r\n",j,nSrcHullPoints);
       //dbug(3,"listLength: %d \r\n",listLength);
       //if ((moduleTypeIndex == 146)||(moduleTypeIndex == 147)) debug_flag=2; else debug_flag = 1;
-      crop_polygon(pixel_R, pixel_U, indexList, &listLength, modules.Height[moduleTypeIndex], modules.Width[moduleTypeIndex], projPoints);
+      crop_polygon(pixel_R, pixel_U, indexList, &listLength, analytic_modules.Height[moduleTypeIndex], analytic_modules.Width[moduleTypeIndex], projPoints);
       //dbug(2,"loopend %d %d\r\n",j,nSrcHullPoints);
       //debug_flag=1;
       //if (listLength>0) dbug(2,"listLength: %d\r\n",listLength);
@@ -1788,7 +1788,7 @@ int any_objects_1(int moduleNumber, void *boundaries)
   int any_objs = 0, i;
   struct height_lims *heightLims = (struct height_lims *) boundaries;
 
-  for(i = 0;i < phantom.numObjects;i++)
+  for(i = 0;i < analytic_phantom.numObjects;i++)
     {
       //dbug(3,"Module number: %d minHeight[%d] = %lf   maxHeight[%d] = %lf     Different?  %d",moduleNumber,i,heightLims.min[i],i,heightLims.max[i],(int)(heightLims.min[i] != heightLims.max[i]));
       if (heightLims[i].min != heightLims[i].max)
@@ -1803,7 +1803,7 @@ int any_objects_2(int moduleNumber, void *boundaries)
   int any_objs = 0, i;
   struct box_lims *boxLims = (struct box_lims *) boundaries;
 
-  for(i = 0;i < phantom.numObjects;i++)
+  for(i = 0;i < analytic_phantom.numObjects;i++)
     {
       if (boxLims[i].minR != boxLims[i].maxR)
 	{any_objs++;break;}
@@ -1817,7 +1817,7 @@ void build_object_list1(double *pix_vlims, int *objectList, int *n_objlist, int
   int i;
   struct height_lims *heightLims = (struct height_lims *) boundaries;
 
-  for(i = 0;i<phantom.numObjects;i++)
+  for(i = 0;i<analytic_phantom.numObjects;i++)
     if((heightLims[i].min <= v+pix_vlims[1]) && (heightLims[i].max >= v+pix_vlims[0]))
       {objectList[n_objlist[0]] = i;n_objlist[0]++;}
 }
@@ -1828,7 +1828,7 @@ void build_object_list2(double *pix_ulims, double *pix_vlims, int *objectList, i
   int i;
   struct box_lims *boxLims = (struct box_lims *) boundaries;
 
-  for(i = 0;i<phantom.numObjects;i++)
+  for(i = 0;i<analytic_phantom.numObjects;i++)
     {
       dbug(3,"if statement parts: %d %d %d %d\n\r",(int) (boxLims[i].minU <= uv[1]+pix_vlims[1]),(int) (boxLims[i].maxU >= uv[1]+pix_vlims[0]),(int) (boxLims[i].minR <= uv[0]+pix_ulims[1]),(int) (boxLims[i].maxR >= uv[0]+pix_ulims[0]));
       dbug(3,"parts for third one: %1.5lf %1.5lf %1.5lf %1.5lf \r\n ",boxLims[0].minR,uv[0],pix_ulims[1],uv[0]+pix_ulims[1]);
@@ -1853,13 +1853,13 @@ void intersections_loop(double *detCenter, double *right, double *up, double *sa
   double *materialBuffer = NULL;
   double *pValueSpectrum = NULL;
 
-  materialBuffer = malloc(sizeof(double)*materials.materialCount);
-  pValueSpectrum = malloc(sizeof(double)*materials.eBinCount);
+  materialBuffer = malloc(sizeof(double)*analytic_materials.materialCount);
+  pValueSpectrum = malloc(sizeof(double)*analytic_materials.eBinCount);
   for(l = 0;l<nSubDets;l++)
     {
       for(i = 0;i<3;i++)
 	subDetCenter[i] = detCenter[i]+right[i]*sampling[2*l]+up[i]*sampling[2*l+1];
-      for(k = 0;k<modules.nSubSources;k++)
+      for(k = 0;k<analytic_modules.nSubSources;k++)
 	{
 	  for(i = 0;i<3;i++)
 	    alpha[i] = subDetCenter[i]-sourcePoints[3*k+i];
@@ -1867,34 +1867,34 @@ void intersections_loop(double *detCenter, double *right, double *up, double *sa
 	  alpha[0] /= rayLength;
 	  alpha[1] /= rayLength;
 	  alpha[2] /= rayLength;
-	  for(i = 0;i<materials.materialCount;i++) materialBuffer[i] = 0.0;
+	  for(i = 0;i<analytic_materials.materialCount;i++) materialBuffer[i] = 0.0;
 	  //dbug(3,"sourcePoints:   [%4.4lf %4.4lf %4.4lf]\r\n",sourcePoints[3*k],sourcePoints[3*k+1],sourcePoints[3*k+2]);
 	  //dbug(3,"alpha:   [%4.4lf %4.4lf %4.4lf]\r\n",alpha[0],alpha[1],alpha[2]);
 	  //dbug(3,"rayLength:   [%4.4lf]\r\n",rayLength);
-	  //dbug(3,"phantom.objectCenter:   [%4.4lf %4.4lf %4.4lf]\r\n",phantom.objectCenter[0],phantom.objectCenter[1],phantom.objectCenter[2]);
+	  //dbug(3,"analytic_phantom.objectCenter:   [%4.4lf %4.4lf %4.4lf]\r\n",analytic_phantom.objectCenter[0],analytic_phantom.objectCenter[1],analytic_phantom.objectCenter[2]);
 	  intersections(&objectList[0],nListObjects,&sourcePoints[3*k],alpha,rayLength,&materialBuffer[0]);
 	  //if (materialBuffer[0]>0.0) dbug(1,"materialBuffer[0]: %4.4lf\r\n",materialBuffer[0]);
-	  for(m = 0;m<materials.eBinCount;m++){
+	  for(m = 0;m<analytic_materials.eBinCount;m++){
 	    pValueSpectrum[m] = 0;
-	    for(n = 0;n<materials.materialCount;n++)
-	      pValueSpectrum[m] += materialBuffer[n]*materials.muTable[n+materials.materialCount*m];
+	    for(n = 0;n<analytic_materials.materialCount;n++)
+	      pValueSpectrum[m] += materialBuffer[n]*analytic_materials.muTable[n+analytic_materials.materialCount*m];
 	  }
 	  //if (pValueSpectrum[0]>0.0) dbug(1,"pValueSpectrum[0]: %4.4lf\r\n",pValueSpectrum[0]);
 
 	  //------------------- Accurate Detector Model, Mingye
 	  if(Accurate_Detector_Model_is_ON)
 	  {
-		  for(m = 0;m < materials.eBinCount;m++)
+		  for(m = 0;m < analytic_materials.eBinCount;m++)
 		  {
-			  int the_index = (detIndex*n_col_oversample + (int)(l/n_row_oversample_add_xtalk))*materials.eBinCount + m;
-			  thisView[the_index] += subviewWeight*detWeights[l]*modules.sourceWeights[k]*exp(-pValueSpectrum[m]);
+			  int the_index = (detIndex*n_col_oversample + (int)(l/n_row_oversample_add_xtalk))*analytic_materials.eBinCount + m;
+			  thisView[the_index] += subviewWeight*detWeights[l]*analytic_modules.sourceWeights[k]*exp(-pValueSpectrum[m]);
 		  }
 	  }
 	  //-------------------
 	  else
 	  {
-		  for(m = 0;m < materials.eBinCount;m++)
-			thisView[detIndex*materials.eBinCount+m] += subviewWeight*detWeights[l]*modules.sourceWeights[k]*exp(-pValueSpectrum[m]);
+		  for(m = 0;m < analytic_materials.eBinCount;m++)
+			thisView[detIndex*analytic_materials.eBinCount+m] += subviewWeight*detWeights[l]*analytic_modules.sourceWeights[k]*exp(-pValueSpectrum[m]);
 	  }
 	}
     }
@@ -1915,24 +1915,24 @@ void set_Accurate_Detector_Model(double *Paras)
 	}
 
 	if(Paras[5]==1 && Paras[1]<3 && Paras[2]==1) //col_xtalk>0, first subview in the first two views 
-		if(modules.Sub[0]>n_col_oversample*n_row_oversample_add_xtalk)
+		if(analytic_modules.Sub[0]>n_col_oversample*n_row_oversample_add_xtalk)
 		{
 			int n1;
 			
-			modules.maxSubDets-=2*n_row_oversample_add_xtalk;
-			modules.Sub[0]-=2*n_row_oversample_add_xtalk;
+			analytic_modules.maxSubDets-=2*n_row_oversample_add_xtalk;
+			analytic_modules.Sub[0]-=2*n_row_oversample_add_xtalk;
 
 			double sum_weights=0;
-			for(n1=0;n1<modules.maxSubDets;n1++)
+			for(n1=0;n1<analytic_modules.maxSubDets;n1++)
 			{
-				modules.Weight[n1]=modules.Weight[n1+n_row_oversample_add_xtalk];
-				sum_weights+=modules.Weight[n1];
+				analytic_modules.Weight[n1]=analytic_modules.Weight[n1+n_row_oversample_add_xtalk];
+				sum_weights+=analytic_modules.Weight[n1];
 			}
-			for(n1=0;n1<modules.maxSubDets;n1++)
-				modules.Weight[n1]/=sum_weights;
+			for(n1=0;n1<analytic_modules.maxSubDets;n1++)
+				analytic_modules.Weight[n1]/=sum_weights;
 
-			for(n1=0;n1<2*modules.maxSubDets;n1++)
-				modules.Sampling[n1]=modules.Sampling[n1+2*n_row_oversample_add_xtalk];
+			for(n1=0;n1<2*analytic_modules.maxSubDets;n1++)
+				analytic_modules.Sampling[n1]=analytic_modules.Sampling[n1+2*n_row_oversample_add_xtalk];
 		}
 
 }
@@ -1954,8 +1954,8 @@ DLLEXPORT void Projector(double *Paras, double subviewWeight, double *thisView,
 	//		Paras[0],Accurate_Detector_Model_is_ON,
 	//		n_col_oversample,n_row_oversample,
 	//		n_col_oversample_add_xtalk,n_row_oversample_add_xtalk,
-	//		modules.Sub[0],modules.maxSubDets,
-	//		modules.Weight[0],modules.Sampling[0]);
+	//		analytic_modules.Sub[0],analytic_modules.maxSubDets,
+	//		analytic_modules.Weight[0],analytic_modules.Sampling[0]);
 	//	fclose(fp);
 	//}
 	//-------------------
@@ -1974,36 +1974,36 @@ DLLEXPORT void Projector(double *Paras, double subviewWeight, double *thisView,
     elapsedTime = (double)( (clk + 0.0) / CLOCKS_PER_SEC);
   }
 
-  objectList = malloc(phantom.numObjects*sizeof(int));
+  objectList = malloc(analytic_phantom.numObjects*sizeof(int));
 
-  switch(modules.moduleOverlapType)
+  switch(analytic_modules.moduleOverlapType)
     {
     case 1:
-      boundaries = malloc(sizeof(struct height_lims)*phantom.numObjects);
+      boundaries = malloc(sizeof(struct height_lims)*analytic_phantom.numObjects);
       break;
     case 2:
-      boundaries = malloc(sizeof(struct box_lims)*phantom.numObjects);
+      boundaries = malloc(sizeof(struct box_lims)*analytic_phantom.numObjects);
       break;
     case 3:
-      //boundaries = malloc(sizeof(struct polygon)*phantom.numObjects);
+      //boundaries = malloc(sizeof(struct polygon)*analytic_phantom.numObjects);
       break;
     }
 
-  //dbug(2,"Number of modules: %d     modules.moduleOverlapType : %d\r\n",nModulesIn, modules.moduleOverlapType);
+  //dbug(2,"Number of analytic_modules: %d     analytic_modules.moduleOverlapType : %d\r\n",nModulesIn, analytic_modules.moduleOverlapType);
   for(moduleNumber = 0;moduleNumber<nModulesIn;moduleNumber++)
     {
       moduleTypeIndex = modTypeInds[moduleNumber];
-      detWeights = &modules.Weight[moduleTypeIndex*modules.maxSubDets];//Weight should be transposed before being passed in
-      UV = &modules.Coords[2*modules.maxPixPerModule*moduleTypeIndex];
+      detWeights = &analytic_modules.Weight[moduleTypeIndex*analytic_modules.maxSubDets];//Weight should be transposed before being passed in
+      UV = &analytic_modules.Coords[2*analytic_modules.maxPixPerModule*moduleTypeIndex];
       //if ((moduleNumber%20)==0) 
-      //  dbug(1,"moduleTypeIndex: %d, modules.Pix[moduleTypeIndex]: %d\n\r",moduleTypeIndex,modules.Pix[moduleTypeIndex]);
-      sampling = &modules.Sampling[2*modules.maxSubDets*moduleTypeIndex];
+      //  dbug(1,"moduleTypeIndex: %d, analytic_modules.Pix[moduleTypeIndex]: %d\n\r",moduleTypeIndex,analytic_modules.Pix[moduleTypeIndex]);
+      sampling = &analytic_modules.Sampling[2*analytic_modules.maxSubDets*moduleTypeIndex];
       //dbug(2,"Module %d  mti:%d\n\r",moduleNumber,moduleTypeIndex);
-      nSubDets = modules.Sub[moduleTypeIndex];
+      nSubDets = analytic_modules.Sub[moduleTypeIndex];
       //dbug(2,"Module %db %d\n\r",moduleNumber,moduleTypeIndex);
       pix_vlims[0] = VERY_BIG;pix_vlims[1] = -VERY_BIG;
       pix_ulims[0] = VERY_BIG;pix_ulims[1] = -VERY_BIG;
-      for(i = 0;i<modules.Sub[moduleTypeIndex];i++){
+      for(i = 0;i<analytic_modules.Sub[moduleTypeIndex];i++){
 	if(sampling[2*i+1]<pix_vlims[0]) pix_vlims[0] = sampling[2*i+1];
 	if(sampling[2*i+1]>pix_vlims[1]) pix_vlims[1] = sampling[2*i+1];
 	if(sampling[2*i]<pix_ulims[0]) pix_ulims[0] = sampling[2*i];
@@ -2017,7 +2017,7 @@ DLLEXPORT void Projector(double *Paras, double subviewWeight, double *thisView,
 
       compute_object_projections(srcHullPoints, nSrcHullPoints, moduleTypeIndex, up, right, center, boundaries);
       //dbug(2,"Module %dc\n\r",moduleNumber);
-      switch(modules.moduleOverlapType)
+      switch(analytic_modules.moduleOverlapType)
 	{
 	case 1:
 	  any_objs = any_objects_1(moduleNumber,boundaries);
@@ -2029,9 +2029,9 @@ DLLEXPORT void Projector(double *Paras, double subviewWeight, double *thisView,
 	  printf("\nERROR: Unrecognized moduleOverlapType!\n\r");
 	}
       if (any_objs) {
-	dbug(1,"There are objects that may project onto module %d.    (%d)\n\r",moduleNumber,modules.moduleOverlapType);
+	dbug(1,"There are objects that may project onto module %d.    (%d)\n\r",moduleNumber,analytic_modules.moduleOverlapType);
         //        dbug(2,"max : %1.12lf   min : %1.12lf\r\n",boundaries[0].max,boundaries[0].min);
-	for(k = 0;k<modules.Pix[moduleTypeIndex];k++){
+	for(k = 0;k<analytic_modules.Pix[moduleTypeIndex];k++){
 	  
 	  // Compute pixel center locations
 	  detIndex = firstDetIndex[moduleNumber]+k;
@@ -2040,12 +2040,12 @@ DLLEXPORT void Projector(double *Paras, double subviewWeight, double *thisView,
 	  
 	  // Build object list
 	  nListObjects = 0;
-	  switch(modules.moduleOverlapType)
+	  switch(analytic_modules.moduleOverlapType)
 	    {
 	    case 1:
 	      build_object_list1(pix_vlims, objectList, &nListObjects, moduleNumber, UV[k*2+1], boundaries);
-	      //if (nListObjects > phantom.numObjects)
-	      //dbug(2,"nListObjects = %d, phantom.numObjects = %d\r\n", nListObjects, phantom.numObjects);
+	      //if (nListObjects > analytic_phantom.numObjects)
+	      //dbug(2,"nListObjects = %d, analytic_phantom.numObjects = %d\r\n", nListObjects, analytic_phantom.numObjects);
 	      //if (nListObjects > 0)
 	      //dbug(2,"nListObjects = %d, k = %d\r\n", nListObjects, k);
 	      break;
@@ -2078,21 +2078,21 @@ DLLEXPORT void Projector(double *Paras, double subviewWeight, double *thisView,
 		  //------------------- Accurate Detector Model, Mingye
 		  if(Accurate_Detector_Model_is_ON)
 		  {
-			  for(m = 0;m<materials.eBinCount;m++)
+			  for(m = 0;m<analytic_materials.eBinCount;m++)
 				  for(l = 0;l<nSubDets;l++)
 				  {
-					  int the_index=(detIndex*n_col_oversample+(int)(l/n_row_oversample_add_xtalk))*materials.eBinCount+m;
+					  int the_index=(detIndex*n_col_oversample+(int)(l/n_row_oversample_add_xtalk))*analytic_materials.eBinCount+m;
 					  thisView[the_index] += subviewWeight*detWeights[l];
 				  }
 		  }
 		  //-------------------
 		  else
 		  {
-			  for(m = 0;m<materials.eBinCount;m++)
-				  thisView[detIndex*materials.eBinCount+m] += subviewWeight;
+			  for(m = 0;m<analytic_materials.eBinCount;m++)
+				  thisView[detIndex*analytic_materials.eBinCount+m] += subviewWeight;
 		  }
 
-	    // thisView[detIndex*materials.eBinCount+m] = thisView[detIndex*materials.eBinCount+m]+subviewWeight[n]*detWeights[l]*sourceWeights;
+	    // thisView[detIndex*analytic_materials.eBinCount+m] = thisView[detIndex*analytic_materials.eBinCount+m]+subviewWeight[n]*detWeights[l]*sourceWeights;
 	  }   // else
 	}    // Loop over pixels (not subpixels)
       }    // if any_objects
@@ -2100,24 +2100,24 @@ DLLEXPORT void Projector(double *Paras, double subviewWeight, double *thisView,
 	//      sum_det_w = 0;                  
 	//      for(l = 0;l<nSubDets;l++) sum_det_w += detWeights[l];
 	// We are assuming the detector weights sum to one and that the sourceWeights sum to one also.
-	for(k = 0;k<modules.Pix[moduleTypeIndex];k++){
+	for(k = 0;k<analytic_modules.Pix[moduleTypeIndex];k++){
 	  detIndex = firstDetIndex[moduleNumber]+k;
 
 	  //------------------- Accurate Detector Model, Mingye
 	  if(Accurate_Detector_Model_is_ON)
 	  {
-		  for(m = 0;m<materials.eBinCount;m++)
+		  for(m = 0;m<analytic_materials.eBinCount;m++)
 			  for(l = 0;l<nSubDets;l++)
 			  {
-				  int the_index=(detIndex*n_col_oversample+(int)(l/n_row_oversample_add_xtalk))*materials.eBinCount+m;
+				  int the_index=(detIndex*n_col_oversample+(int)(l/n_row_oversample_add_xtalk))*analytic_materials.eBinCount+m;
 				  thisView[the_index] += subviewWeight*detWeights[l];
 			  }
 	  }
 	  //-------------------
 	  else
 	  {
-		  for(m = 0;m<materials.eBinCount;m++)
-			  thisView[detIndex*materials.eBinCount+m] += subviewWeight;
+		  for(m = 0;m<analytic_materials.eBinCount;m++)
+			  thisView[detIndex*analytic_materials.eBinCount+m] += subviewWeight;
 	  }
 
 	}
@@ -2125,10 +2125,10 @@ DLLEXPORT void Projector(double *Paras, double subviewWeight, double *thisView,
       if (debug_flag>0)
 	{
 	  min_val = VERY_BIG;
-	  for(k = 0;k<modules.Pix[moduleTypeIndex];k++){
+	  for(k = 0;k<analytic_modules.Pix[moduleTypeIndex];k++){
 	    detIndex = firstDetIndex[moduleNumber]+k;
-	    for(m = 0;m<materials.eBinCount;m++)
-	      if (thisView[detIndex*materials.eBinCount+m]<min_val) min_val = thisView[detIndex*materials.eBinCount+m];
+	    for(m = 0;m<analytic_materials.eBinCount;m++)
+	      if (thisView[detIndex*analytic_materials.eBinCount+m]<min_val) min_val = thisView[detIndex*analytic_materials.eBinCount+m];
 	  }
 	  dbug(1,"minval: %f (module #%d)\n\r", min_val,moduleNumber);
 	}
@@ -2232,8 +2232,8 @@ DLLEXPORT void Projector_threaded(double *Paras, double subviewWeight, double *t
 	//		Paras[0],Accurate_Detector_Model_is_ON,
 	//		n_col_oversample,n_row_oversample,
 	//		n_col_oversample_add_xtalk,n_row_oversample_add_xtalk,
-	//		modules.Sub[0],modules.maxSubDets,
-	//		modules.Weight[0],modules.Sampling[0]);
+	//		analytic_modules.Sub[0],analytic_modules.maxSubDets,
+	//		analytic_modules.Weight[0],analytic_modules.Sampling[0]);
 	//	fclose(fp);
 	//}
 	//-------------------
diff --git a/gecatsim/clib_build/src/analytic_projector.h b/gecatsim/clib_build/src/analytic_projector.h
index f24a485..aae8b99 100644
--- a/gecatsim/clib_build/src/analytic_projector.h
+++ b/gecatsim/clib_build/src/analytic_projector.h
@@ -27,7 +27,7 @@ int n_row_oversample = 1;
 int n_col_oversample_add_xtalk = 1;
 int n_row_oversample_add_xtalk = 1;
 
-struct module_info
+struct analytic_module_info
 {
   double *Height;
   double *Width;
@@ -43,7 +43,7 @@ struct module_info
   int moduleOverlapType;
 };
 
-struct module_info modules = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0,0,0,0};
+struct analytic_module_info analytic_modules = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0,0,0,0};
 
 struct pnt
 {
@@ -61,7 +61,7 @@ struct bounding_info
 
 struct bounding_info bounding = {NULL,NULL};
 
-struct phantom_info
+struct analytic_phantom_info
 {
   int numObjects;
   int *objectType;
@@ -77,7 +77,7 @@ struct phantom_info
   double xbounds[2];
 };
 
-struct phantom_info phantom = {0,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,{0,0}};
+struct analytic_phantom_info analytic_phantom = {0,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,{0,0}};
 
 struct material_info
 {
@@ -86,7 +86,7 @@ struct material_info
   double *muTable;
 };
 
-struct material_info materials = {0,0,NULL};
+struct material_info analytic_materials = {0,0,NULL};
 
 struct projector_args
 {
diff --git a/gecatsim/clib_build/src/phantom.tab.cpp b/gecatsim/clib_build/src/phantom.tab.cpp
index 56e2034..0d324c0 100644
--- a/gecatsim/clib_build/src/phantom.tab.cpp
+++ b/gecatsim/clib_build/src/phantom.tab.cpp
@@ -64,7 +64,7 @@
 typedef struct phantomlval {
   double val;
   double vec[3];
-  char str[128];
+  char str[2048];
 } PHANTOMLVAL;
 
 #endif
@@ -290,7 +290,7 @@ typedef struct phantomlval {
 typedef union YYSTYPE {
   double val;
   double vec[3];
-  char str[128];
+  char str[2048];
 } YYSTYPE;
 /* Line 191 of yacc.c.  */
 //#line 283 "phantom.tab.cpp"
diff --git a/gecatsim/clib_build/src/phantom.tab.hpp b/gecatsim/clib_build/src/phantom.tab.hpp
index 446bc69..cb368ec 100644
--- a/gecatsim/clib_build/src/phantom.tab.hpp
+++ b/gecatsim/clib_build/src/phantom.tab.hpp
@@ -110,7 +110,7 @@
 typedef union YYSTYPE {
   double val;
   double vec[3];
-  char str[128];
+  char str[2048];
 } YYSTYPE;
 /* Line 1275 of yacc.c.  */
 #line 113 "phantom.tab.hpp"
diff --git a/gecatsim/pyfiles/C_Projector_SetData.py b/gecatsim/pyfiles/C_Projector_SetData.py
index bc8e77f..e5fa8c5 100644
--- a/gecatsim/pyfiles/C_Projector_SetData.py
+++ b/gecatsim/pyfiles/C_Projector_SetData.py
@@ -31,6 +31,8 @@ def C_Projector_SetData(cfg, viewId, subViewId):
 def get_projector_id(cfg):
     if not isinstance(cfg.phantom.callback, list):
         phantom_callbacks = [cfg.phantom.callback]
+    elif isinstance(cfg.phantom.callback, list) and isinstance(cfg.phantom.callback[0], list):
+        phantom_callbacks = cfg.phantom.callback[cfg.viewId]
     else:
         phantom_callbacks = cfg.phantom.callback
     projectorIDs = []
diff --git a/gecatsim/pyfiles/C_Projector_Voxelized.py b/gecatsim/pyfiles/C_Projector_Voxelized.py
index 13ea98b..46f55ed 100644
--- a/gecatsim/pyfiles/C_Projector_Voxelized.py
+++ b/gecatsim/pyfiles/C_Projector_Voxelized.py
@@ -51,7 +51,8 @@ def C_Projector_Voxelized(cfg, viewId, subViewId):
             unused7 = 7
             freeTheMemory = 0
             
-            if viewId==cfg.sim.stopViewId and subViewId==cfg.sim.subViewCount-1 \
+            # ZJY recalcPht is helpful in reducing memory usage in dynamic phantom
+            if (cfg.physics.recalcPht or (viewId==cfg.sim.stopViewId and subViewId==cfg.sim.subViewCount-1))\
                 and srcId==src.nSamples-1 and matId==cfg.phantom.numberOfMaterials-1:
                 freeTheMemory = 1
             
diff --git a/gecatsim/pyfiles/CommonTools.py b/gecatsim/pyfiles/CommonTools.py
index 3f74878..59b6a9d 100644
--- a/gecatsim/pyfiles/CommonTools.py
+++ b/gecatsim/pyfiles/CommonTools.py
@@ -26,11 +26,12 @@ def make_col(a):
     return a
 
 
-def feval(func_name, *args):
+def feval(funcName, *args):
     module_paths = [
-        func_name,
-        f"gecatsim.pyfiles.{func_name}",
-        f"gecatsim.pyfiles.FlatPanel.{func_name}"
+        funcName,
+        f"gecatsim.pyfiles.{funcName}",
+        f"gecatsim.pyfiles.FlatPanel.{funcName}",
+        f"gecatsim.src.app.pyfiles.{funcName}"
     ]
 
     module = None
@@ -42,15 +43,15 @@ def feval(func_name, *args):
             continue
 
     if module is None:
-        raise ImportError(f"Could not import module: {func_name}. Tried paths: {module_paths}")
+        raise ImportError(f"Could not import module: {funcName}. Tried paths: {module_paths}")
 
     # Extract the function name (the part after the last dot)
-    func_name_only = func_name.split('.')[-1]
+    funcName_only = funcName.split('.')[-1]
 
     try:
-        func = getattr(module, func_name_only)
+        func = getattr(module, funcName_only)
     except AttributeError:
-        raise AttributeError(f"Function '{func_name_only}' not found in module '{module.__name__}'")
+        raise AttributeError(f"Function '{funcName_only}' not found in module '{module.__name__}'")
 
     return func(*args)
 
@@ -74,16 +75,35 @@ def load_C_lib():
 
 
 class emptyCFG:
-    pass
+    
+    def __str__(self):
+        attrs = ', '.join(f"{key}={value!r}" for key, value in self.__dict__.items())
+        return f"{self.__class__.__name__}({attrs})"
+    
+    def __repr__(self):
+        return self.__str__()
+    
+
 
+import sys
 
 class PathHelper:
     def __init__(self):
         self._base_dir = os.getcwd()
         # Locate paths of lib and data.
         self.paths = {}
-        self.paths["main"] = os.path.dirname(os.path.abspath(__file__))
-        self.paths["top"] = os.path.split(self.paths["main"])[0]
+        
+        # Handle PyInstaller bundled environment
+        if getattr(sys, 'frozen', False) and hasattr(sys, '_MEIPASS'):
+            # Running as compiled executable
+            bundle_dir = sys._MEIPASS
+            self.paths["main"] = os.path.join(bundle_dir, 'app', 'pyfiles')
+            self.paths["top"] = os.path.join(bundle_dir, 'app')
+        else:
+            # Running as normal Python script
+            self.paths["main"] = os.path.dirname(os.path.abspath(__file__))
+            self.paths["top"] = os.path.split(self.paths["main"])[0]
+        
         self.paths["bowtie"] = os.path.join(self.paths["top"], 'bowtie')
         self.paths["cfg"] = os.path.join(self.paths["top"], 'cfg')
         self.paths["dose"] = os.path.join(self.paths["top"], 'dose')
@@ -259,7 +279,6 @@ def source_cfg(*para):
         
     # initialize structs in cfg and structs
     attrList = ['sim', 'det', 'detNew', 'src', 'srcNew', 'spec', 'protocol', 'scanner', 'phantom', 'physics', 'recon', 'dose']
-    
     for attr in attrList:
         # Ensure cfg has the attribute
         if not hasattr(cfg, attr):
@@ -381,27 +400,60 @@ def overlap2d(oldimg, old_pos_x, old_pos_y, pos_x, pos_y):
 def rawread(fname, dataShape, dataType):
     # dataType is for numpy, ONLY allows: 'float'/'single', 'double', 'int'/'int32', 'uint'/'uint32', 'int8', 'int16' 
     #          they are single, double, int32, uin32, int8, int16
-    with open(fname, 'rb') as fin:
-        data = fin.read()
-    
     # https://docs.python.org/3/library/struct.html
-    switcher = {'float': ['f', 4, np.single], 
-                'single': ['f', 4, np.single], 
-                'double': ['d', 8, np.double], 
+    switcher = {'float': ['f', 4, np.float32], 
+                'float32': ['f', 4, np.float32],
+                'single': ['f', 4, np.float32], 
+                'double': ['d', 8, np.float64], 
                 'int': ['i', 4, np.int32], 
                 'uint': ['I', 4, np.uint32],  
                 'int32': ['i', 4, np.int32], 
                 'uint32': ['I', 4, np.uint32], 
+                'uint8': ['B', 1, np.uint8],
                 'int8': ['b', 1, np.int8], 
                 'int16': ['h', 2, np.int16]}
-    fmt = switcher[dataType]
-    data = struct.unpack("%d%s" % (len(data)/fmt[1], fmt[0]), data)
     
-    data = np.array(data, dtype=fmt[2])
+    fmt_char, byte_size, np_dtype = switcher[dataType]
+    
+    # Use numpy.fromfile for direct binary reading - much faster
+    try:
+        data = np.fromfile(fname, dtype=np_dtype)
+    except:
+        # Fallback to manual reading if fromfile fails
+        with open(fname, 'rb') as fin:
+            raw_data = fin.read()
+        num_elements = len(raw_data) // byte_size
+        data = struct.unpack(f"{num_elements}{fmt_char}", raw_data)
+        data = np.array(data, dtype=np_dtype)
+    
     if dataShape:
         data = data.reshape(dataShape)
     
     return data
+    
+    # with open(fname, 'rb') as fin:
+    #     data = fin.read()
+    
+    # # https://docs.python.org/3/library/struct.html
+    # switcher = {'float': ['f', 4, np.single], 
+    #             'float32': ['f', 4, np.float32],
+    #             'single': ['f', 4, np.single], 
+    #             'double': ['d', 8, np.double], 
+    #             'int': ['i', 4, np.int32], 
+    #             'uint': ['I', 4, np.uint32],  
+    #             'int32': ['i', 4, np.int32], 
+    #             'uint32': ['I', 4, np.uint32], 
+    #             'uint8': ['B', 1, np.uint8],
+    #             'int8': ['b', 1, np.int8], 
+    #             'int16': ['h', 2, np.int16]}
+    # fmt = switcher[dataType]
+    # data = struct.unpack("%d%s" % (len(data)/fmt[1], fmt[0]), data)
+    
+    # data = np.array(data, dtype=fmt[2])
+    # if dataShape:
+    #     data = data.reshape(dataShape)
+    
+    # return data
 
 def rawwrite(fname, data):
     with open(fname, 'wb') as fout:
diff --git a/gecatsim/pyfiles/Detector_RayAngles_2D.py b/gecatsim/pyfiles/Detector_RayAngles_2D.py
index 562024c..e3d15eb 100644
--- a/gecatsim/pyfiles/Detector_RayAngles_2D.py
+++ b/gecatsim/pyfiles/Detector_RayAngles_2D.py
@@ -30,7 +30,7 @@ def Detector_RayAngles_2D(cfg):
 
         xyzDet = nm.repmat(det.modCoords[m, :], nCells, 1) + \
             det.cellCoords @ (np.c_[det.uvecs[m, :], det.vvecs[m, :]].T)
-        xyzR = xyzDet - nm.repmat(xyzSrc, nCells, 1);
+        xyzR = xyzDet - nm.repmat(xyzSrc, nCells, 1)
 
         cellInd = range(startInd, startInd+nCells)
         rayDistance[cellInd] = vectornorm(xyzR.T)
diff --git a/gecatsim/pyfiles/PhantomProjectorWrapper.py b/gecatsim/pyfiles/PhantomProjectorWrapper.py
index 1a912f8..ed519fb 100644
--- a/gecatsim/pyfiles/PhantomProjectorWrapper.py
+++ b/gecatsim/pyfiles/PhantomProjectorWrapper.py
@@ -1,28 +1,168 @@
+gecatsim\pyfiles\PhantomProjectorWrapper.py
+@@ -1,466 +1,94 @@
+"""
+PhantomProjectorWrapper.py - Phantom and Projector Configuration Management
+
+=============================================================================
+OVERVIEW - WHERE THIS FITS IN THE XCIST PIPELINE
+=============================================================================
+
+This module is called during the phantom_scan phase of CatSim.run_all():
+
+    CatSim.run_all()
+      └── phantom_scan()
+            └── one_scan()
+                  └── RunModels() with callback_map including:
+                        - 'Phantom': PhantomWrapper()     [initialization phase]
+                        - 'Projector': ProjectorWrapper() [per-view projection]
+
+The module handles THREE phantom configuration scenarios:
+
+1. SINGLE PHANTOM (simplest case):
+   cfg.phantom.projectorCallback = "C_Projector_Voxelized"
+   → One phantom, same for all views
+
+2. MULTIPLE PHANTOMS PER VIEW (e.g., breast + lesion + calcifications):
+   cfg.phantom.projectorCallback = ["C_Projector_Voxelized", "Analytic_Projector", ...]
+   cfg.phantom.filename = ["breast.json", "lesion.json", ...]
+   → Multiple phantoms combined, same set for all views
+
+3. DYNAMIC PHANTOM (nested list - different phantom per view):
+   cfg.phantom.projectorCallback = [
+       ["C_Projector_Voxelized"],           # view 0
+       ["C_Projector_Voxelized"],           # view 1
+       ["Analytic_Projector", "C_Proj..."], # view 2 (multiple phantoms)
+       ...
+   ]
+   → Each view can have different phantom configuration (e.g., motion simulation)
+
+=============================================================================
+MEMORY CONSIDERATIONS
+=============================================================================
+
+For multiple phantoms per view, partial_subview is allocated:
+    Shape: [nSamples, totalNumCells, nEbin]
+    Size: nSamples × totalNumCells × nEbin × 4 bytes
+    Example: 9 × 170,572,500 × 74 × 4 = ~450 GB for full mammo resolution!
+
+This is the LARGEST memory allocation in the entire pipeline and occurs
+because each phantom's transmission must be computed independently before
+combining (multiplicative combination for overlapping phantoms).
+
+For single phantom, no partial_subview is needed - projection writes
+directly to cfg.thisSubView.
+
+=============================================================================
+FUNCTION CALL SEQUENCE
+=============================================================================
+
+Per-view execution flow:
+
+    RunModels(cfg, viewId, subViewId)
+        │
+        ├── [First view only] PhantomWrapper(cfg)
+        │       └── Handles dynamic phantom: extracts cfg.phantom for viewId
+        │
+        └── ProjectorWrapper(cfg, viewId, subViewId)
+                │
+                ├── [Dynamic phantom] Re-extract cfg.phantom for viewId
+                │
+                ├── [Multiple phantoms] 
+                │   ├── SplitCfgPhantom() → list of individual phantom cfgs
+                │   ├── Allocate partial_subview [nSamples, nCells, nEbin]
+                │   ├── Loop over each phantom:
+                │   │   ├── CopyCfgPhantom() → copy phantom config
+                │   │   ├── feval(phantom.callback) → Phantom_Voxelized()
+                │   │   └── feval(projectorCallback) → C_Projector_Voxelized()
+                │   └── Combine: thisSubView *= sum(partial_subview)
+                │
+                └── [Single phantom]
+                    ├── feval(phantom.callback) → Phantom_Voxelized()
+                    └── feval(projectorCallback) → C_Projector_Voxelized()
+
+=============================================================================
+"""
+
 import sys
 from copy import copy, deepcopy
 from gecatsim.pyfiles.CommonTools import *
 
 def CopyCfgPhantom(cfgfrom, cfgto=None):
+    """
+    Copy phantom configuration from one CFG object to another.
+    
+    Used during multiple-phantom iteration to update the main cfg with
+    each individual phantom's configuration while preserving all other
+    cfg attributes (detector, source, spectrum, etc.).
+    
+    Args:
+        cfgfrom: Source CFG object containing phantom to copy
+        cfgto: Destination CFG object (created if None)
+    
+    Returns:
+        cfgto with phantom attributes deep-copied from cfgfrom
+    
+    Note: Uses deepcopy to avoid reference issues when phantom attributes
+    are modified during projection.
+    """
     if cfgto is None:
         cfgto = CFG()
     cfgto.phantom = deepcopy(cfgfrom.phantom)
 
     return cfgto
 
-# split cfg.phantom from a list into 
 def SplitCfgPhantom(cfg):
+    """
+    Split a multi-phantom configuration into a list of single-phantom CFGs.
+    
+    When cfg.phantom contains lists (multiple phantoms to combine), this
+    function creates separate CFG objects for each phantom so they can be
+    processed individually by their respective projectors.
+    
+    Example input:
+        cfg.phantom.filename = ["breast.raw", "lesion.raw", "calc.raw"]
+        cfg.phantom.projectorCallback = ["C_Projector", "C_Projector", "Analytic"]
+        cfg.phantom.materialList = [["adipose"], ["tumor"], ["calcium"]]
+    
+    Output:
+        cfg_list[0].phantom.filename = "breast.raw"
+        cfg_list[0].phantom.projectorCallback = "C_Projector"
+        cfg_list[0].phantom.materialList = ["adipose"]
+        ...and so on for each phantom
+    
+    Args:
+        cfg: Main CFG with phantom attributes as lists
+    
+    Returns:
+        List of CFG objects, each with single phantom configuration
+    
+    Note: Only phantom attributes are copied; other cfg attributes (det, src, etc.)
+    are accessed from the main cfg via CopyCfgPhantom during iteration.
+    """
     cfg_list = []
     cfgphantom_attrs = [x for x in dir(cfg.phantom) if not x.startswith('__')]
     for i in range(len(cfg.phantom.filename)):
         thiscfg = CFG()
         for thisattr in cfgphantom_attrs:
-            setattr(thiscfg.phantom, thisattr, deepcopy(getattr(cfg.phantom, thisattr)[i]))
+            # broadcast in this case
+            if not isinstance(getattr(cfg.phantom, thisattr), list) or len(getattr(cfg.phantom, thisattr))==1:
+                setattr(thiscfg.phantom, thisattr, deepcopy(getattr(cfg.phantom, thisattr)))
+            else:
+                setattr(thiscfg.phantom, thisattr, deepcopy(getattr(cfg.phantom, thisattr)[i]))
 
         cfg_list.append(thiscfg)
 
     return cfg_list
 
 def PhantomWrapper(cfg):
+    # nested loop, like [[pro1, proj2], [proj1], [proj2], ...]
+    is_dynamic_phantom = False
+    if isinstance(cfg.phantom.callback, list) and isinstance(cfg.phantom.callback[0], list):
+        is_dynamic_phantom = True
+        _cfg_list = SplitCfgPhantom(cfg)
+        _origCfgPhantom = CopyCfgPhantom(cfg)
+        cfg = CopyCfgPhantom(_cfg_list[cfg.viewId], cfg)
+
     if isinstance(cfg.phantom.callback, list):
         # if duplicate phantom callback, raise an error
         if (len(cfg.phantom.callback) != len(set(cfg.phantom.callback))) and len(cfg.phantom.callback) > 1:
@@ -33,6 +173,7 @@ def PhantomWrapper(cfg):
         origCfgPhantom = CopyCfgPhantom(cfg)
         cfg_nom_list = []
         for thiscfg in cfg_list:
+            # use cfg here because we want to use cfg's other attr than phantom without the heavy copy overhead
             cfg = CopyCfgPhantom(thiscfg, cfg)
             cfg = feval(cfg.phantom.callback, cfg)
             if hasattr(cfg.phantom, 'numberOfMaterials'):
@@ -43,14 +184,27 @@ def PhantomWrapper(cfg):
                 cfg_nom_list.append(None)
 
         # restore original CFG
-        cfg = CopyCfgPhantom(origCfgPhantom, cfg)
-        cfg.phantom.numberOfMaterials = cfg_nom_list
+        if is_dynamic_phantom:
+            cfg = CopyCfgPhantom(_origCfgPhantom, cfg)
+            cfg.phantom.numberOfMaterials = [None]*len(cfg.phantom.callback)
+            cfg.phantom.numberOfMaterials[cfg.viewId] = cfg_nom_list
+        else:
+            cfg = CopyCfgPhantom(origCfgPhantom, cfg)
+            cfg.phantom.numberOfMaterials = cfg_nom_list
     else:
         cfg = feval(cfg.phantom.callback, cfg)
 
+
     return cfg
 
 def ProjectorWrapper(cfg, viewId, subViewId):
+    is_dynamic_phantom = False
+    if isinstance(cfg.phantom.callback, list) and isinstance(cfg.phantom.callback[0], list):
+        is_dynamic_phantom = True
+        _cfg_list = SplitCfgPhantom(cfg)
+        _origCfgPhantom = CopyCfgPhantom(cfg)
+        cfg = CopyCfgPhantom(_cfg_list[cfg.viewId], cfg)
+
     if isinstance(cfg.phantom.projectorCallback, list):
         cfg_list = SplitCfgPhantom(cfg)
         origCfgPhantom = CopyCfgPhantom(cfg)
@@ -62,4 +216,7 @@ def ProjectorWrapper(cfg, viewId, subViewId):
     else:
         cfg = feval(cfg.phantom.projectorCallback, cfg, viewId, subViewId)
 
+    if is_dynamic_phantom:
+        cfg = CopyCfgPhantom(_origCfgPhantom, cfg)
+
     return cfg
diff --git a/gecatsim/pyfiles/Phantom_Analytic.py b/gecatsim/pyfiles/Phantom_Analytic.py
index fc7aeb4..f5350b2 100644
--- a/gecatsim/pyfiles/Phantom_Analytic.py
+++ b/gecatsim/pyfiles/Phantom_Analytic.py
@@ -648,10 +648,12 @@ def Phantom_Analytic_SetObjects(objs=None, debug=False):
             phantObject[i]['s'] = objs['clip'][i][:,3].T
         else:
             #TODO: consider change to []
-            #phantObject[i]['eta'] = []
-            #phantObject[i]['s'] = []
-            phantObject[i]['eta'] = None
-            phantObject[i]['s'] = None
+            # phantObject[i]['eta'] = []
+            # phantObject[i]['s'] = []
+            # phantObject[i]['eta'] = None
+            # phantObject[i]['s'] = None
+            phantObject[i]['eta'] = np.empty((3, 0))
+            phantObject[i]['s'] = np.empty(0)
         p1=pars[6]*np.pi / 180
         p2=pars[7]*np.pi / 180
         p3=pars[8]*np.pi / 180
diff --git a/gecatsim/pyfiles/Xray_Filter.py b/gecatsim/pyfiles/Xray_Filter.py
index dadf7b5..a50dbd3 100644
--- a/gecatsim/pyfiles/Xray_Filter.py
+++ b/gecatsim/pyfiles/Xray_Filter.py
@@ -20,49 +20,59 @@ def Xray_Filter(cfg):
 
 def flat_filter(cfg):
     '''
-    Apply the transmittance of flat filter at source side, dim: [Ebin, totalNumCells]
+    Apply the transmittance of flat filter at source side, optimized version
     '''
-    cosineFactors = 1/np.cos(cfg.det.gammas)/np.cos(cfg.det.alphas)
+    cosineFactors = 1 / (np.cos(cfg.det.gammas) * np.cos(cfg.det.alphas))
     
     Evec = cfg.sim.Evec
-    trans = np.ones([cfg.det.totalNumCells, cfg.spec.nEbin], dtype = np.single)
+    trans = np.ones([cfg.det.totalNumCells, cfg.spec.nEbin], dtype=np.single)
     
     if hasattr(cfg.protocol, "flatFilter"):
-        for ii in range(0, round(len(cfg.protocol.flatFilter)/2)):
+        # Pre-allocate array for mu values
+        mu_total = np.zeros_like(Evec, dtype=np.single)
+        
+        for ii in range(0, len(cfg.protocol.flatFilter) // 2):
             material = cfg.protocol.flatFilter[2*ii]
             depth = cfg.protocol.flatFilter[2*ii+1]
             mu = GetMu(material, Evec)
-            trans *= np.exp(-depth*0.1*cosineFactors @ mu.reshape(1, mu.size))
-    cfg.src.filterTrans *= trans
+            mu_total += depth * 0.1 * mu
+        
+        # Single exponential calculation instead of multiple
+        trans *= np.exp(-np.outer(cosineFactors, mu_total))
     
+    cfg.src.filterTrans *= trans
     return cfg
 
 def bowtie_filter(cfg):
     '''
-    Apply the transmittance of bowtie filter, dim: [Ebin, totalNumCells]
+    Apply the transmittance of bowtie filter, optimized version
     '''
     if not cfg.protocol.bowtie:
         return cfg
     
-    # find bowtie file
     bowtieFile = my_path.find("bowtie", cfg.protocol.bowtie, ".txt")
-
-    # read bowtie file
     data = np.loadtxt(bowtieFile, dtype=np.single, comments=['#', '%'])
     
     gammas0 = data[:, 0]
-    t0 = data[:, 1:] # thickness in cm
+    t0 = data[:, 1:]
     bowtieMaterials = ['Al', 'graphite', 'Cu', 'Ti']
     
     Evec = cfg.spec.Evec
     gammas1 = cfg.det.gammas
+    cos_alphas_inv = 1.0 / np.cos(cfg.det.alphas)
+    
+    # Pre-compute all mu values
+    mu_array = np.array([GetMu(material, Evec) for material in bowtieMaterials], dtype=np.single)
+    
+    # Vectorized interpolation for all materials at once
+    f_interps = [interpolate.interp1d(gammas0, t0[:, i], kind='linear', fill_value='extrapolate', assume_sorted=True) 
+                 for i in range(len(bowtieMaterials))]
     
-    muT = 0
-    for i in range(len(bowtieMaterials)):
-        mu = GetMu(bowtieMaterials[i], Evec)
-        f = interpolate.interp1d(gammas0, t0[:, i], kind='linear', fill_value='extrapolate')
-        t1 = f(gammas1)/np.cos(cfg.det.alphas)
-        muT += t1 @ mu.reshape(1, mu.size)
+    # Compute total attenuation in one operation
+    muT = np.zeros((len(gammas1), len(Evec)), dtype=np.single)
+    for i, f in enumerate(f_interps):
+        t1 = f(gammas1) * cos_alphas_inv
+        muT += np.outer(t1, mu_array[i])
     
     trans = np.exp(-muT)
     cfg.src.filterTrans *= trans