diff --git a/cleedpy/cleed/src/leed.c b/cleedpy/cleed/src/leed.c
index 23c3dbf..ca95350 100644
--- a/cleedpy/cleed/src/leed.c
+++ b/cleedpy/cleed/src/leed.c
@@ -1,5 +1,28 @@
 #include "leed.h"
 
+void print_phase_shift(struct phs_str phs_shift)
+{
+    int i;
+    printf("Phase shift:\n");
+    printf("  lmax: %d\n", phs_shift.lmax);
+    printf("  neng: %d\n", phs_shift.neng);
+    printf("  t_type: %d\n", phs_shift.t_type);
+    printf("  eng_max: %lf\n", phs_shift.eng_max);
+    printf("  eng_min: %lf\n", phs_shift.eng_min);
+    printf("  energy: %lf\n", phs_shift.energy[0]);
+    printf("  dr: %lf\n", phs_shift.dr[0]);
+    printf("  input_file: %s\n", phs_shift.input_file);
+    printf("  pshift: ");
+    for (i=0; i<phs_shift.neng; i++)
+    {
+        printf("%lf ", phs_shift.pshift[i]);
+        if (i % 10 == 0)
+            printf("\n");
+    }
+    printf("\n");
+
+}
+
 real ** leed(
     char * par_file,
     char * bul_file,
@@ -18,7 +41,7 @@ real ** leed(
 
     int n_beams_now, n_beams_set;
 
-    int i_c, i_set, offset;
+    int i_c, i_set, offset, i;
     int i_layer;
     int energy_index;
     int n_set;
@@ -40,11 +63,12 @@ real ** leed(
     FILE *res_stream;
 
     // Read input parameters
-    printf("HERE %s\n", bul_file);
     inp_rdbul_nd(&bulk, &phs_shifts, bul_file);
     inp_rdpar(&v_par, &eng, bulk, bul_file);
     inp_rdovl_nd(&over, &phs_shifts, bulk, par_file);
     inp_showbop(bulk, over, phs_shifts);
+    for (i=0; (phs_shifts + i)->lmax != I_END_OF_LIST; i++)
+        print_phase_shift(phs_shifts[i]);
 
 
     // Construct energy list
@@ -55,9 +79,6 @@ real ** leed(
         energy_list[energy_index] = eng->ini + energy_index * eng->stp;
     }
 
-    // Printing stuff
-    inp_showbop(bulk, over, phs_shifts);
-
 
     res_stream = fopen(res_file,"w");
 
diff --git a/cleedpy/cleed/src/libmat/matalloc.c b/cleedpy/cleed/src/libmat/matalloc.c
index 51e2fd0..3cb9bc3 100644
--- a/cleedpy/cleed/src/libmat/matalloc.c
+++ b/cleedpy/cleed/src/libmat/matalloc.c
@@ -54,232 +54,161 @@ mat matalloc(mat M, int rows, int cols, int num_type)
 
 *********************************************************************/
 {
-int mat_type;
-size_t no_of_elts;
-real *ptr, *ptr_end;
+    int mat_type;
+    size_t no_of_elts;
+    real *ptr, *ptr_end;
 
-#ifdef CONTROL
- fprintf(STDCTR,"(matalloc): enter function rows: %d cols: %d num_type: %d\n",
-         rows, cols, num_type);
-#endif
-
-/*********************************************************************
-  Check the validity of the pointer M and cols/rows first
-*********************************************************************/
-
- if(matcheck(M) < 0)
- {
-#ifdef ERROR
-   fprintf(STDERR," *** error (matalloc): Invalid pointer \n");
-#endif
-#ifdef EXIT_ON_ERROR
-   exit(1);
-#else
-   return(NULL);
-#endif
- }
+    /*********************************************************************
+     Check the validity of the pointer M and cols/rows first
+    *********************************************************************/
 
- if( (cols < 1) || (rows < 1) )
- {
-#ifdef ERROR
-   fprintf(STDERR," *** error (matalloc): Invalid number of cols/rows \n");
-#endif
-#ifdef EXIT_ON_ERROR
-   exit(1);
-#else
-   return(NULL);
-#endif
- }
-
-/*********************************************************************
-  Find the correct matrix type (high byte of num_type) and number
-  types (low byte of num_type).
-*********************************************************************/
+    if(matcheck(M) < 0)
+    {
+        fprintf(STDERR," *** error (matalloc): Invalid pointer \n");
+        exit(1);
+    }
 
- mat_type = num_type & MAT_MASK;
-
- if(mat_type == 0)
- {
-   if(cols != rows) mat_type = MAT_NORMAL;
-   else
-     if(cols != 1)  mat_type = MAT_SQUARE;
-     else           mat_type = MAT_SCALAR;
- }
- else if( (mat_type == MAT_DIAG) && (cols != rows) )
- {
-#ifdef ERROR
-   fprintf(STDERR," *** error (matalloc): different numbers of cols ");
-   fprintf(STDERR,"and rows are incompatible with diagonal matrix type\n");
-#endif
-#ifdef EXIT_ON_ERROR
-   exit(1);
-#else
-   return(NULL);
-#endif
- }
+    if( (cols < 1) || (rows < 1) )
+    {
+        fprintf(STDERR," *** error (matalloc): Invalid number of cols/rows \n");
+        exit(1);
+    }
 
- num_type = num_type & NUM_MASK;
+    /*********************************************************************
+     Find the correct matrix type (high byte of num_type) and number
+    types (low byte of num_type).
+    *********************************************************************/
 
-#ifdef CONTROL_X
- fprintf(STDCTR,"(matalloc): num_type: %x mat_type: %x \n", num_type, mat_type);
-#endif
+    mat_type = num_type & MAT_MASK;
 
-/*********************************************************************
-  If M points to the right matrix type already, only reset all matrix
-  elements.
-*********************************************************************/
+    if(mat_type == 0)
+    {
+        if(cols != rows) mat_type = MAT_NORMAL;
+        else
+            if(cols != 1) mat_type = MAT_SQUARE;
+            else mat_type = MAT_SCALAR;
+    }
+    else if( (mat_type == MAT_DIAG) && (cols != rows) )
+    {
+        fprintf(STDERR," *** error (matalloc): different numbers of cols and rows are incompatible with diagonal matrix type\n");
+        exit(1);
+    }
 
- if( (matcheck(M) > 0) &&
-     (M->cols == cols) &&
-     (M->rows == rows) &&
-     (M->num_type == num_type) &&
-     (M->mat_type == mat_type) )
- {
-   if(M->num_type == NUM_COMPLEX)
-   {
-   /*
-     M->mat_type < MAT_DIAG means square, normal or scalar
-   */
-     if(M->mat_type < MAT_DIAG)
-     {
-       ptr_end = M->rel + M->rows*M->cols;
-       for(ptr = M->rel; ptr <= ptr_end; ptr ++) *ptr = 0.;
-
-       ptr_end = M->iel + M->rows*M->cols;
-       for(ptr = M->iel; ptr <= ptr_end; ptr ++) *ptr = 0.;
-     }
-     if(M->mat_type == MAT_DIAG)
-     {
-       ptr_end = M->rel + M->cols;
-       for(ptr = M->rel; ptr <= ptr_end; ptr ++) *ptr = 0.;
-
-       ptr_end = M->iel + M->cols;
-       for(ptr = M->iel; ptr <= ptr_end; ptr ++) *ptr = 0.;
-     }
-   } /* NUM_COMPLEX */
-
-   if(M->num_type == NUM_REAL)
-   {
-     if(M->mat_type < MAT_DIAG)
-     {
-       ptr_end = M->rel + M->rows*M->cols;
-       for(ptr = M->rel; ptr <= ptr_end; ptr ++) *ptr = 0.;
-     }
-     if(M->mat_type == MAT_DIAG)
-     {
-       ptr_end = M->rel + M->cols;
-       for(ptr = M->rel; ptr <= ptr_end; ptr ++) *ptr = 0.;
-     }
-   } /* NUM_REAL */
- /*
-   Set magic number and return
- */
-   M->mag_no = MATRIX;
-   return(M);
- }  /* reset existing matrix */
+    num_type = num_type & NUM_MASK;
+    /*******************************************************************************
+     If M points to the right matrix type already, only reset all matrix elements.
+    *******************************************************************************/
 
-/*********************************************************************
-  If not: allocate new memory
+    if( (matcheck(M) > 0) && (M->cols == cols) && (M->rows == rows) && (M->num_type == num_type) && (M->mat_type == mat_type) )
+    {
+        if(M->num_type == NUM_COMPLEX)
+        {
+            // M->mat_type < MAT_DIAG means square, normal or scalar
+            if(M->mat_type < MAT_DIAG)
+            {
+                ptr_end = M->rel + M->rows*M->cols;
+                for(ptr = M->rel; ptr <= ptr_end; ptr ++) *ptr = 0.;
+
+                ptr_end = M->iel + M->rows*M->cols;
+                for(ptr = M->iel; ptr <= ptr_end; ptr ++) *ptr = 0.;
+            }
+            if(M->mat_type == MAT_DIAG)
+            {
+                ptr_end = M->rel + M->cols;
+                for(ptr = M->rel; ptr <= ptr_end; ptr ++) *ptr = 0.;
+
+                ptr_end = M->iel + M->cols;
+                for(ptr = M->iel; ptr <= ptr_end; ptr ++) *ptr = 0.;
+            }
+        } /* NUM_COMPLEX */
+
+        if(M->num_type == NUM_REAL)
+        {
+            if(M->mat_type < MAT_DIAG)
+            {
+            ptr_end = M->rel + M->rows*M->cols;
+            for(ptr = M->rel; ptr <= ptr_end; ptr ++) *ptr = 0.;
+            }
+            if(M->mat_type == MAT_DIAG)
+            {
+            ptr_end = M->rel + M->cols;
+            for(ptr = M->rel; ptr <= ptr_end; ptr ++) *ptr = 0.;
+            }
+        } /* NUM_REAL */
+        // Set magic number and return
+        M->mag_no = MATRIX;
+        return(M);
+    }  /* reset existing matrix */
+
+    /*********************************************************************
+    If not: allocate new memory
+
+    - if M == NULL: allocate a new matrix header, blk_type = BLK_SINGLE
+    - if M != NULL: update the old header, use former blk_type.
+    *********************************************************************/
+
+    if( M == NULL)
+    {
+        M = ( mat )malloc(sizeof(struct mat_str));
+        M->blk_type = BLK_SINGLE;
+    }
+    else  // M != NULL
+    {
+        if (M->iel != NULL) free(M->iel);
+        if (M->rel != NULL) free(M->rel);
+    }
 
-  - if M == NULL: allocate a new matrix header, blk_type = BLK_SINGLE
-  - if M != NULL: update the old header, use former blk_type.
-*********************************************************************/
+    M->cols = cols;
+    M->rows = rows;
+    M->num_type = num_type;
+    M->mat_type = mat_type;
 
- if( M == NULL)
- {
-#ifdef CONTROL_X
-   fprintf(STDCTR,"(matalloc): create new matrix structure\n");
-#endif
-   M = ( mat )malloc(sizeof(struct mat_str));
-   M->blk_type = BLK_SINGLE;
- }
- else  /* M != NULL */
- {
-#ifdef CONTROL_X
-   fprintf(STDCTR,"(matalloc): reuse old matrix structure\n");
-#endif
-   if (M->iel != NULL) free(M->iel);
-   if (M->rel != NULL) free(M->rel);
- }
-
- M->cols = cols;
- M->rows = rows;
- M->num_type = num_type;
- M->mat_type = mat_type;
-
- if (M->mat_type == MAT_DIAG) no_of_elts = M->cols+1;
- else                         no_of_elts = M->rows * M->cols + 1;
-
- switch(num_type)
- {
-  case(NUM_REAL):
-  {
-#ifdef CONTROL_X
-    fprintf(STDCTR,"(matalloc): allocate %d real matrix elements\n", no_of_elts);
-#endif
-    M->iel = NULL;
-    M->rel = (real*)calloc( no_of_elts, sizeof(real));
+    if(M->mat_type == MAT_DIAG)
+        no_of_elts = M->cols+1;
+    else
+        no_of_elts = M->rows * M->cols + 1;
 
-    if (M->rel == NULL)
+    switch(num_type)
     {
-      free(M);
-#ifdef ERROR
-     fprintf(STDERR,"*** error (matalloc) allocation error\n");
-#endif
-#ifdef EXIT_ON_ERROR
-     exit(1);
-#else
-     return(NULL);
-#endif
-    }
-    break;
-  }  /* NUM_REAL */
-
-  case(NUM_COMPLEX):
-  {
-#ifdef CONTROL_X
-    fprintf(STDCTR,"(matalloc): allocate 2 * %d complex matrix elements\n",
-            no_of_elts);
-#endif
-    M->rel = (real*)calloc( no_of_elts, sizeof(real));
-    M->iel = (real*)calloc( no_of_elts, sizeof(real));
+        case(NUM_REAL):
+        {
+            M->iel = NULL;
+            M->rel = (real*)calloc( no_of_elts, sizeof(real));
+            if (M->rel == NULL)
+            {
+                free(M);
+                fprintf(STDERR,"*** error (matalloc) allocation error\n");
+                exit(1);
+            }
+            break;
+        }  /* NUM_REAL */
+
+        case(NUM_COMPLEX):
+        {
+            M->rel = (real*)calloc( no_of_elts, sizeof(real));
+            M->iel = (real*)calloc( no_of_elts, sizeof(real));
+
+            if( (M->rel == NULL) || (M->iel == NULL) )
+            {
+                if (M->rel != NULL) free(M->rel);
+                if (M->iel != NULL) free(M->iel);
+                free(M);
+                fprintf(STDERR,"*** error (matalloc) allocation error\n");
+                exit(1);
+            }
+            break;
+        }  /* NUM_COMPLEX */
+
+        default: // Unknown number type.
+        {
+            fprintf(STDERR, " *** error (matalloc) wrong number type: %d\n", num_type);
+            exit(1);
+
+        }
+    } /* switch */
+
+    // Set magic number and return
+    M->mag_no = MATRIX;
+    return(M);
 
-    if( (M->rel == NULL) || (M->iel == NULL) )
-    {
-      if (M->rel != NULL) free(M->rel);
-      if (M->iel != NULL) free(M->iel);
-      free(M);
-#ifdef ERROR
-      fprintf(STDERR,"*** error (matalloc) allocation error\n");
-#endif
-#ifdef EXIT_ON_ERROR
-      exit(1);
-#else
-      return(NULL);
-#endif
-    }
-    break;
-  }  /* NUM_COMPLEX */
-
-  default:   /* unknown number type */
-  {
-#ifdef ERROR
-    fprintf(STDERR,
-            " *** error (matalloc) wrong number type: %d\n", num_type);
-#endif
-#ifdef EXIT_ON_ERROR
-     exit(1);
-#else
-     return(NULL);
-#endif
-  }
- } /* switch */
-/*
-  Set magic number and return
-*/
-#ifdef CONTROL_X
-    fprintf(STDCTR,"(matalloc): set magic number and return\n\n");
-#endif
- M->mag_no = MATRIX;
- return(M);
 }  /* end of function matalloc */
diff --git a/cleedpy/cleed/src/linpphasend.c b/cleedpy/cleed/src/linpphasend.c
index e871b3a..1ee884c 100644
--- a/cleedpy/cleed/src/linpphasend.c
+++ b/cleedpy/cleed/src/linpphasend.c
@@ -265,8 +265,8 @@ real faux;
  phs_shifts->lmax = lmax;
  nl = lmax + 1;
 
- phs_shifts->energy = (real *)malloc( neng * sizeof(real) );
- phs_shifts->pshift = (real *)malloc( neng * nl * sizeof(real) );
+ phs_shifts->energy = (real *)calloc( neng, sizeof(real) );
+ phs_shifts->pshift = (real *)calloc( neng * nl, sizeof(real) );
 
 
  for( i_eng = 0;