From 68d776b5e78bbe7f3b757b3b9e31e6243ebb4fc3 Mon Sep 17 00:00:00 2001
From: Julian Kauth <julian.kauth@fau.de>
Date: Sun, 17 Jul 2022 13:50:13 +0200
Subject: [PATCH] remove copy from phi_in to psi_out in evolve_wave

Instead of copying phi_in to psi_out value by value in evolve_wave_half we
can redesign how we use our pointers and instead of copying everything we
can simply swap pointers afterwards. In some tests this resulted in a speed
up of ~10%
---
 mangrove.c        | 26 +++++++++-----------------
 wave_3d.c         | 26 +++++++++-----------------
 wave_billiard.c   | 37 ++++++++++++++++++++-----------------
 wave_comparison.c | 29 +++++++++--------------------
 wave_energy.c     | 26 +++++++++-----------------
 5 files changed, 56 insertions(+), 88 deletions(-)

diff --git a/mangrove.c b/mangrove.c
index e27c8f6..32235f8 100644
--- a/mangrove.c
+++ b/mangrove.c
@@ -453,7 +453,7 @@ void init_bc_phase(double left_bc[NY], double top_bc[NX], double bot_bc[NX])
 // //     printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]);
 // }
 
-void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], double *psi_out[NX], 
+void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX],
                       short int *xy_in[NX])
 /* time step of field evolution */
 /* phi is value of field at time t, psi at time t-1 */
@@ -512,7 +512,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
 
                 /* evolve phi */
                 phi_out[i][j] = -y + 2*x + tcc[i][j]*delta - KAPPA*x - tgamma[i][j]*(x-y);
-                psi_out[i][j] = x;
             }
         }
     }
@@ -557,7 +556,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[0][j] = x;
         }
     }
     
@@ -593,7 +591,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[NX-1][j] = x;
         }
     }
     
@@ -656,7 +653,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[i][NY-1] = x;
         }
     }
     
@@ -719,7 +715,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[i][0] = x;
         }
     }
     
@@ -737,20 +732,19 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
             {
                 if (phi_out[i][j] > VMAX) phi_out[i][j] = VMAX;
                 if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX;
-                if (psi_out[i][j] > VMAX) psi_out[i][j] = VMAX;
-                if (psi_out[i][j] < -VMAX) psi_out[i][j] = -VMAX;
             }
         }
     }
 }
 
 
-void evolve_wave(double *phi[NX], double *psi[NX], double *phi_tmp[NX], double *psi_tmp[NX], short int *xy_in[NX])
+void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX])
 /* time step of field evolution */
 /* phi is value of field at time t, psi at time t-1 */
 {
-    evolve_wave_half(phi, psi, phi_tmp, psi_tmp, xy_in);
-    evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in);
+    evolve_wave_half(phi, psi, tmp, xy_in);
+    evolve_wave_half(tmp, phi, psi, xy_in);
+    evolve_wave_half(psi, tmp, phi, xy_in);
 }
 
 
@@ -848,7 +842,7 @@ void animation()
 {
     double time, scale, diss, rgb[3], hue, y, dissip, ej, gradient[2], dx, dy, dt, xleft, xright, 
             length, fx, fy, force[2];
-    double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX];
+    double *phi[NX], *psi[NX], *tmp[NX];
     short int *xy_in[NX], redraw = 0;
     int i, j, k, n, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q;
     static int imin, imax;
@@ -862,8 +856,7 @@ void animation()
     {
         phi[i] = (double *)malloc(NY*sizeof(double));
         psi[i] = (double *)malloc(NY*sizeof(double));
-        phi_tmp[i] = (double *)malloc(NY*sizeof(double));
-        psi_tmp[i] = (double *)malloc(NY*sizeof(double));
+        tmp[i] = (double *)malloc(NY*sizeof(double));
         xy_in[i] = (short int *)malloc(NY*sizeof(short int));
     }
     mangrove = (t_mangrove *)malloc(NMAXCIRCLES*sizeof(t_mangrove));    /* mangroves */  
@@ -977,7 +970,7 @@ void animation()
         for (j=0; j<NVID; j++) 
         {
 //             printf("%i ", j);
-            evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
+            evolve_wave(phi, psi, tmp, xy_in);
 //             if (i % 10 == 9) oscillate_linear_wave(0.2*scale, 0.15*(double)(i*NVID + j), -1.5, YMIN, -1.5, YMAX, phi, psi);
         }
         
@@ -1258,8 +1251,7 @@ void animation()
     {
         free(phi[i]);
         free(psi[i]);
-        free(phi_tmp[i]);
-        free(psi_tmp[i]);
+        free(tmp[i]);
         free(xy_in[i]);
     }
     free(mangrove);
diff --git a/wave_3d.c b/wave_3d.c
index e283aff..64591e3 100644
--- a/wave_3d.c
+++ b/wave_3d.c
@@ -266,7 +266,7 @@ FILE *time_series_left, *time_series_right;
 double courant2, courantb2;  /* Courant parameters squared */
 
 
-void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out[NX*NY], double psi_out[NX*NY], 
+void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out[NX*NY],
                       short int xy_in[NX*NY], double tc[NX*NY], double tcc[NX*NY], double tgamma[NX*NY])
 // void evolve_wave_half(double *phi_in, double *psi_in, double *phi_out, double *psi_out, 
 //                       short int *xy_in[NX])
@@ -295,7 +295,6 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
 
                 /* evolve phi */
                 phi_out[i*NY+j] = -y + 2*x + tcc[i*NY+j]*delta - KAPPA*x - tgamma[i*NY+j]*(x-y);
-                psi_out[i*NY+j] = x;
             }
         }
     }
@@ -333,7 +332,6 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
                     break;
                 }
             }
-            psi_out[j] = x;
         }
     }
     
@@ -369,7 +367,6 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
                     break;
                 }
             }
-            psi_out[(NX-1)*NY+j] = x;
         }
     }
     
@@ -419,7 +416,6 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
                     break;
                 }
             }
-            psi_out[i*NY+NY-1] = x;
         }
     }
     
@@ -469,7 +465,6 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
                     break;
                 }
             }
-            psi_out[i*NY] = x;
         }
     }
     
@@ -487,21 +482,20 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
             {
                 if (phi_out[i*NY+j] > VMAX) phi_out[i*NY+j] = VMAX;
                 if (phi_out[i*NY+j] < -VMAX) phi_out[i*NY+j] = -VMAX;
-                if (psi_out[i*NY+j] > VMAX) psi_out[i*NY+j] = VMAX;
-                if (psi_out[i*NY+j] < -VMAX) psi_out[i*NY+j] = -VMAX;
             }
         }
     }
 }
 
 
-void evolve_wave(double phi[NX*NY], double psi[NX*NY], double phi_tmp[NX*NY], double psi_tmp[NX*NY], short int xy_in[NX*NY],
+void evolve_wave(double phi[NX*NY], double psi[NX*NY], double tmp[NX*NY], short int xy_in[NX*NY],
     double tc[NX*NY], double tcc[NX*NY], double tgamma[NX*NY])
 /* time step of field evolution */
 /* phi is value of field at time t, psi at time t-1 */
 {
-    evolve_wave_half(phi, psi, phi_tmp, psi_tmp, xy_in, tc, tcc, tgamma);
-    evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in, tc, tcc, tgamma);
+    evolve_wave_half(phi, psi, tmp, xy_in, tc, tcc, tgamma);
+    evolve_wave_half(tmp, phi, psi, xy_in, tc, tcc, tgamma);
+    evolve_wave_half(psi, tmp, phi, xy_in, tc, tcc, tgamma);
 }
 
 
@@ -542,7 +536,7 @@ void viewpoint_schedule(int i)
 void animation()
 {
     double time, scale, ratio, startleft[2], startright[2], sign, r2, xy[2], fade_value; 
-    double *phi, *psi, *phi_tmp, *psi_tmp, *total_energy, *color_scale, *tc, *tcc, *tgamma;
+    double *phi, *psi, *phi_tmp, *tmp, *total_energy, *color_scale, *tc, *tcc, *tgamma;
     short int *xy_in;
     int i, j, s, sample_left[2], sample_right[2], period = 0, fade;
     static int counter = 0;
@@ -559,8 +553,7 @@ void animation()
     xy_in = (short int *)malloc(NX*NY*sizeof(short int));
     phi = (double *)malloc(NX*NY*sizeof(double));
     psi = (double *)malloc(NX*NY*sizeof(double));
-    phi_tmp = (double *)malloc(NX*NY*sizeof(double));
-    psi_tmp = (double *)malloc(NX*NY*sizeof(double));
+    tmp = (double *)malloc(NX*NY*sizeof(double));
     total_energy = (double *)malloc(NX*NY*sizeof(double));
     color_scale = (double *)malloc(NX*NY*sizeof(double));
     tc = (double *)malloc(NX*NY*sizeof(double));
@@ -663,7 +656,7 @@ void animation()
         draw_wave_3d(0, phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
         for (j=0; j<NVID; j++) 
         {
-            evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in, tc, tcc, tgamma);
+            evolve_wave(phi, psi, tmp, xy_in, tc, tcc, tgamma);
             if (SAVE_TIME_SERIES)
             {
                 wave_value = (long int)(phi[sample_left[0]*NY+sample_left[1]]*1.0e16);
@@ -777,8 +770,7 @@ void animation()
     free(xy_in);
     free(phi);
     free(psi);
-    free(phi_tmp);
-    free(psi_tmp);
+    free(tmp);
     free(total_energy);
     free(color_scale);
     free(tc);
diff --git a/wave_billiard.c b/wave_billiard.c
index 97508d5..307ad94 100644
--- a/wave_billiard.c
+++ b/wave_billiard.c
@@ -227,7 +227,7 @@ double courant2, courantb2;  /* Courant parameters squared */
 
 
 
-void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], double *psi_out[NX], 
+void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX],
                       short int *xy_in[NX])
 /* time step of field evolution */
 /* phi is value of field at time t, psi at time t-1 */
@@ -277,7 +277,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
 
                 /* evolve phi */
                 phi_out[i][j] = -y + 2*x + tcc[i][j]*delta - KAPPA*x - tgamma[i][j]*(x-y);
-                psi_out[i][j] = x;
             }
         }
     }
@@ -315,7 +314,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[0][j] = x;
         }
     }
     
@@ -351,7 +349,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[NX-1][j] = x;
         }
     }
     
@@ -401,7 +398,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[i][NY-1] = x;
         }
     }
     
@@ -451,7 +447,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[i][0] = x;
         }
     }
     
@@ -469,20 +464,30 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
             {
                 if (phi_out[i][j] > VMAX) phi_out[i][j] = VMAX;
                 if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX;
-                if (psi_out[i][j] > VMAX) psi_out[i][j] = VMAX;
-                if (psi_out[i][j] < -VMAX) psi_out[i][j] = -VMAX;
             }
         }
     }
 }
 
 
-void evolve_wave(double *phi[NX], double *psi[NX], double *phi_tmp[NX], double *psi_tmp[NX], short int *xy_in[NX])
+void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX])
 /* time step of field evolution */
 /* phi is value of field at time t, psi at time t-1 */
 {
-    evolve_wave_half(phi, psi, phi_tmp, psi_tmp, xy_in);
-    evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in);
+    // For the purpose of these comments w[t], w[t-1], w[t+1] are used to refer
+    // to phi, psi and the result respectively to avoid confusion with the
+    // passed parameter names.
+    // At the beginning w[t] is saved in phi, w[t-1] in psi and tmp is space
+    // for the next wave state w[t+1]. Take w[t] and w[t-1] to calculate the
+    // next wave state. Write this new state in temp
+    evolve_wave_half(phi, psi, tmp, xy_in);
+    // now w[t] is saved in tmp, w[t-1] in phi and the result is written to psi
+    evolve_wave_half(tmp, phi, psi, xy_in);
+    // now w[t] is saved in psi, w[t-1] in tmp and the result is written to phi
+    evolve_wave_half(psi, tmp, phi, xy_in);
+    // now w[t] is saved in phi, w[t-1] in psi and tmp is free again to take
+    // the new wave state w[t+1] in the next call to this function, thus
+    // matching the given parameter names again
 }
 
 
@@ -514,7 +519,7 @@ void draw_color_bar_palette(int plot, double range, int palette, int fade, doubl
 void animation()
 {
     double time, scale, ratio, startleft[2], startright[2], sign, r2, xy[2], fade_value; 
-    double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX], *total_energy[NX], *color_scale[NX];
+    double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *color_scale[NX];
     short int *xy_in[NX];
     int i, j, s, sample_left[2], sample_right[2], period = 0, fade;
     static int counter = 0;
@@ -531,8 +536,7 @@ void animation()
     {
         phi[i] = (double *)malloc(NY*sizeof(double));
         psi[i] = (double *)malloc(NY*sizeof(double));
-        phi_tmp[i] = (double *)malloc(NY*sizeof(double));
-        psi_tmp[i] = (double *)malloc(NY*sizeof(double));
+        tmp[i] = (double *)malloc(NY*sizeof(double));
         total_energy[i] = (double *)malloc(NY*sizeof(double));
         xy_in[i] = (short int *)malloc(NY*sizeof(short int));
         color_scale[i] = (double *)malloc(NY*sizeof(double));
@@ -654,7 +658,7 @@ void animation()
         else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, i, PLOT, COLOR_PALETTE, 0, 1.0);
         for (j=0; j<NVID; j++) 
         {
-            evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
+            evolve_wave(phi, psi, tmp, xy_in);
             if (SAVE_TIME_SERIES)
             {
                 wave_value = (long int)(phi[sample_left[0]][sample_left[1]]*1.0e16);
@@ -760,8 +764,7 @@ void animation()
     {
         free(phi[i]);
         free(psi[i]);
-        free(phi_tmp[i]);
-        free(psi_tmp[i]);
+        free(tmp[i]);
         free(total_energy[i]);
         free(xy_in[i]);
         free(color_scale[i]);
diff --git a/wave_comparison.c b/wave_comparison.c
index 2cb7ef7..d3f627a 100644
--- a/wave_comparison.c
+++ b/wave_comparison.c
@@ -232,7 +232,7 @@ double courant2, courantb2;  /* Courant parameters squared */
 /* animation part    */
 /*********************/
 
-void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], double *psi_out[NX], 
+void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX],
                       short int *xy_in[NX])
 /* time step of field evolution */
 /* phi is value of field at time t, psi at time t-1 */
@@ -280,7 +280,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
 
                 /* evolve phi */
                 phi_out[i][j] = -y + 2*x + tcc[i][j]*delta - KAPPA*x - tgamma[i][j]*(x-y);
-                psi_out[i][j] = x;
             }
         }
         for (j=jmid+1; j<NY-1; j++){
@@ -293,7 +292,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
 
                 /* evolve phi */
                 phi_out[i][j] = -y + 2*x + tcc[i][j]*delta - KAPPA*x - tgamma[i][j]*(x-y);
-                psi_out[i][j] = x;
             }
         }
     }
@@ -340,7 +338,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[0][j] = x;
         }
     }
     
@@ -382,7 +379,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[NX-1][j] = x;
         }
     }
     
@@ -440,7 +436,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[i][jmid-1] = x;
         }
     }
     
@@ -498,7 +493,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[i][0] = x;
         }
     }
     
@@ -556,7 +550,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[i][NY-1] = x;
         }
     }
     
@@ -614,7 +607,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[i][jmid] = x;
         }
     }
     
@@ -634,8 +626,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
             {
                 if (phi_out[i][j] > VMAX) phi_out[i][j] = VMAX;
                 if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX;
-                if (psi_out[i][j] > VMAX) psi_out[i][j] = VMAX;
-                if (psi_out[i][j] < -VMAX) psi_out[i][j] = -VMAX;
             }
         }
     }
@@ -643,12 +633,13 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
 }
 
 
-void evolve_wave(double *phi[NX], double *psi[NX], double *phi_tmp[NX], double *psi_tmp[NX], short int *xy_in[NX])
+void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX])
 /* time step of field evolution */
 /* phi is value of field at time t, psi at time t-1 */
 {
-    evolve_wave_half(phi, psi, phi_tmp, psi_tmp, xy_in);
-    evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in);
+    evolve_wave_half(phi, psi, tmp, xy_in);
+    evolve_wave_half(tmp, phi, psi, xy_in);
+    evolve_wave_half(psi, tmp, phi, xy_in);
 }
 
 
@@ -663,7 +654,7 @@ void draw_color_bar(int plot, double range)
 void animation()
 {
     double time, scale, energies[6], top_energy, bottom_energy;
-    double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX];
+    double *phi[NX], *psi[NX], *tmp[NX];
     short int *xy_in[NX];
     int i, j, s, counter = 0;
 
@@ -672,8 +663,7 @@ void animation()
     {
         phi[i] = (double *)malloc(NY*sizeof(double));
         psi[i] = (double *)malloc(NY*sizeof(double));
-        phi_tmp[i] = (double *)malloc(NY*sizeof(double));
-        psi_tmp[i] = (double *)malloc(NY*sizeof(double));
+        tmp[i] = (double *)malloc(NY*sizeof(double));
         xy_in[i] = (short int *)malloc(NY*sizeof(short int));
     }
     
@@ -763,7 +753,7 @@ void animation()
         
         for (j=0; j<NVID; j++) 
         {
-            evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
+            evolve_wave(phi, psi, tmp, xy_in);
 //             if (i % 10 == 9) oscillate_linear_wave(0.2*scale, 0.15*(double)(i*NVID + j), -1.5, YMIN, -1.5, YMAX, phi, psi);
         }
         
@@ -833,8 +823,7 @@ void animation()
     {
         free(phi[i]);
         free(psi[i]);
-        free(phi_tmp[i]);
-        free(psi_tmp[i]);
+        free(tmp[i]);
         free(xy_in[i]);
     }
 
diff --git a/wave_energy.c b/wave_energy.c
index 03c3e28..c23868a 100644
--- a/wave_energy.c
+++ b/wave_energy.c
@@ -515,7 +515,7 @@ void evolve_wave_half_old(double *phi_in[NX], double *psi_in[NX], double *phi_ou
 }
 
 
-void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], double *psi_out[NX], 
+void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX],
                       short int *xy_in[NX])
 /* time step of field evolution */
 /* phi is value of field at time t, psi at time t-1 */
@@ -563,7 +563,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
 
                 /* evolve phi */
                 phi_out[i][j] = -y + 2*x + tcc[i][j]*delta - KAPPA*x - tgamma[i][j]*(x-y);
-                psi_out[i][j] = x;
             }
         }
     }
@@ -601,7 +600,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[0][j] = x;
         }
     }
     
@@ -637,7 +635,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[NX-1][j] = x;
         }
     }
     
@@ -688,7 +685,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                    break;
                 }
             }
-            psi_out[i][jmid-1] = x;
         }
     }
     
@@ -737,7 +733,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                     break;
                 }
             }
-            psi_out[i][0] = x;
         }
     }
     
@@ -755,8 +750,6 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
             {
                 if (phi_out[i][j] > VMAX) phi_out[i][j] = VMAX;
                 if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX;
-                if (psi_out[i][j] > VMAX) psi_out[i][j] = VMAX;
-                if (psi_out[i][j] < -VMAX) psi_out[i][j] = -VMAX;
             }
         }
     }
@@ -764,14 +757,15 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
 }
 
 
-void evolve_wave(double *phi[NX], double *psi[NX], double *phi_tmp[NX], double *psi_tmp[NX], short int *xy_in[NX])
+void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX])
 /* time step of field evolution */
 /* phi is value of field at time t, psi at time t-1 */
 {
 //     evolve_wave_half_old(phi, psi, phi_tmp, psi_tmp, xy_in);
 //     evolve_wave_half_old(phi_tmp, psi_tmp, phi, psi, xy_in);
-    evolve_wave_half(phi, psi, phi_tmp, psi_tmp, xy_in);
-    evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in);
+    evolve_wave_half(phi, psi, tmp, xy_in);
+    evolve_wave_half(tmp, phi, psi, xy_in);
+    evolve_wave_half(psi, tmp, phi, xy_in);
 }
 
 
@@ -779,7 +773,7 @@ void evolve_wave(double *phi[NX], double *psi[NX], double *phi_tmp[NX], double *
 void animation()
 {
     double time, scale, energies[6], top_energy, bottom_energy;
-    double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX];
+    double *phi[NX], *psi[NX], *tmp[NX];
     short int *xy_in[NX];
     int i, j, s;
 
@@ -788,8 +782,7 @@ void animation()
     {
         phi[i] = (double *)malloc(NY*sizeof(double));
         psi[i] = (double *)malloc(NY*sizeof(double));
-        phi_tmp[i] = (double *)malloc(NY*sizeof(double));
-        psi_tmp[i] = (double *)malloc(NY*sizeof(double));
+        tmp[i] = (double *)malloc(NY*sizeof(double));
         xy_in[i] = (short int *)malloc(NY*sizeof(short int));
     }
     
@@ -852,7 +845,7 @@ void animation()
         
         for (j=0; j<NVID; j++) 
         {
-            evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
+            evolve_wave(phi, psi, tmp, xy_in);
 //             if (i % 10 == 9) oscillate_linear_wave(0.2*scale, 0.15*(double)(i*NVID + j), -1.5, YMIN, -1.5, YMAX, phi, psi);
         }
         
@@ -884,8 +877,7 @@ void animation()
     {
         free(phi[i]);
         free(psi[i]);
-        free(phi_tmp[i]);
-        free(psi_tmp[i]);
+        free(tmp[i]);
         free(xy_in[i]);
     }