daviddoji/YouTube-simulations
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #5 from JulianKauth/main

Browse Source 
Remove copy from phi_in to psi_out in evolve_wave_half
This commit is contained in:
Nils Berglund 2022-07-19 19:07:27 +02:00 committed by GitHub
commit db761031f2
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 56 additions and 88 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
mangrove.c
View File

@ -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]); // // 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]) short int *xy_in[NX])
/* time step of field evolution */ /* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */ /* 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 */ /* evolve phi */
phi_out[i][j] = -y + 2*x + tcc[i][j]*delta - KAPPA*x - tgamma[i][j]*(x-y); 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; 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; 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; 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; 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 (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 */ /* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */ /* 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, psi, tmp, xy_in);
evolve_wave_half(phi_tmp, psi_tmp, phi, psi, 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, double time, scale, diss, rgb[3], hue, y, dissip, ej, gradient[2], dx, dy, dt, xleft, xright,
length, fx, fy, force[2]; 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; short int *xy_in[NX], redraw = 0;
int i, j, k, n, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q; int i, j, k, n, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q;
static int imin, imax; static int imin, imax;
@ -862,8 +856,7 @@ void animation()
{ {
phi[i] = (double *)malloc(NY*sizeof(double)); phi[i] = (double *)malloc(NY*sizeof(double));
psi[i] = (double *)malloc(NY*sizeof(double)); psi[i] = (double *)malloc(NY*sizeof(double));
phi_tmp[i] = (double *)malloc(NY*sizeof(double)); tmp[i] = (double *)malloc(NY*sizeof(double));
psi_tmp[i] = (double *)malloc(NY*sizeof(double));
xy_in[i] = (short int *)malloc(NY*sizeof(short int)); xy_in[i] = (short int *)malloc(NY*sizeof(short int));
} }
mangrove = (t_mangrove *)malloc(NMAXCIRCLES*sizeof(t_mangrove)); /* mangroves */ mangrove = (t_mangrove *)malloc(NMAXCIRCLES*sizeof(t_mangrove)); /* mangroves */
@ -977,7 +970,7 @@ void animation()
for (j=0; j<NVID; j++) for (j=0; j<NVID; j++)
{ {
// printf("%i ", 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); // 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(phi[i]);
free(psi[i]); free(psi[i]);
free(phi_tmp[i]); free(tmp[i]);
free(psi_tmp[i]);
free(xy_in[i]); free(xy_in[i]);
} }
free(mangrove); free(mangrove);

26
wave_3d.c
View File

@ -266,7 +266,7 @@ FILE *time_series_left, *time_series_right;
double courant2, courantb2; /* Courant parameters squared */ 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]) 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, // void evolve_wave_half(double *phi_in, double *psi_in, double *phi_out, double *psi_out,
// short int *xy_in[NX]) // 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 */ /* evolve phi */
phi_out[i*NY+j] = -y + 2*x + tcc[i*NY+j]*delta - KAPPA*x - tgamma[i*NY+j]*(x-y); 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; 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; 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; 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; 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 (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]) double tc[NX*NY], double tcc[NX*NY], double tgamma[NX*NY])
/* time step of field evolution */ /* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */ /* 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, psi, tmp, xy_in, tc, tcc, tgamma);
evolve_wave_half(phi_tmp, psi_tmp, phi, psi, 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() void animation()
{ {
double time, scale, ratio, startleft[2], startright[2], sign, r2, xy[2], fade_value; 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; short int *xy_in;
int i, j, s, sample_left[2], sample_right[2], period = 0, fade; int i, j, s, sample_left[2], sample_right[2], period = 0, fade;
static int counter = 0; static int counter = 0;
@ -559,8 +553,7 @@ void animation()
xy_in = (short int *)malloc(NX*NY*sizeof(short int)); xy_in = (short int *)malloc(NX*NY*sizeof(short int));
phi = (double *)malloc(NX*NY*sizeof(double)); phi = (double *)malloc(NX*NY*sizeof(double));
psi = (double *)malloc(NX*NY*sizeof(double)); psi = (double *)malloc(NX*NY*sizeof(double));
phi_tmp = (double *)malloc(NX*NY*sizeof(double)); tmp = (double *)malloc(NX*NY*sizeof(double));
psi_tmp = (double *)malloc(NX*NY*sizeof(double));
total_energy = (double *)malloc(NX*NY*sizeof(double)); total_energy = (double *)malloc(NX*NY*sizeof(double));
color_scale = (double *)malloc(NX*NY*sizeof(double)); color_scale = (double *)malloc(NX*NY*sizeof(double));
tc = (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); draw_wave_3d(0, phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
for (j=0; j<NVID; j++) 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) if (SAVE_TIME_SERIES)
{ {
wave_value = (long int)(phi[sample_left[0]*NY+sample_left[1]]*1.0e16); wave_value = (long int)(phi[sample_left[0]*NY+sample_left[1]]*1.0e16);
@ -777,8 +770,7 @@ void animation()
free(xy_in); free(xy_in);
free(phi); free(phi);
free(psi); free(psi);
free(phi_tmp); free(tmp);
free(psi_tmp);
free(total_energy); free(total_energy);
free(color_scale); free(color_scale);
free(tc); free(tc);

37
wave_billiard.c
View File

@ -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]) short int *xy_in[NX])
/* time step of field evolution */ /* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */ /* 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 */ /* evolve phi */
phi_out[i][j] = -y + 2*x + tcc[i][j]*delta - KAPPA*x - tgamma[i][j]*(x-y); 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; 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; 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; 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; 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 (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 */ /* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */ /* phi is value of field at time t, psi at time t-1 */
{ {
evolve_wave_half(phi, psi, phi_tmp, psi_tmp, xy_in); // For the purpose of these comments w[t], w[t-1], w[t+1] are used to refer
evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in); // 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() void animation()
{ {
double time, scale, ratio, startleft[2], startright[2], sign, r2, xy[2], fade_value; 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]; short int *xy_in[NX];
int i, j, s, sample_left[2], sample_right[2], period = 0, fade; int i, j, s, sample_left[2], sample_right[2], period = 0, fade;
static int counter = 0; static int counter = 0;
@ -531,8 +536,7 @@ void animation()
{ {
phi[i] = (double *)malloc(NY*sizeof(double)); phi[i] = (double *)malloc(NY*sizeof(double));
psi[i] = (double *)malloc(NY*sizeof(double)); psi[i] = (double *)malloc(NY*sizeof(double));
phi_tmp[i] = (double *)malloc(NY*sizeof(double)); tmp[i] = (double *)malloc(NY*sizeof(double));
psi_tmp[i] = (double *)malloc(NY*sizeof(double));
total_energy[i] = (double *)malloc(NY*sizeof(double)); total_energy[i] = (double *)malloc(NY*sizeof(double));
xy_in[i] = (short int *)malloc(NY*sizeof(short int)); xy_in[i] = (short int *)malloc(NY*sizeof(short int));
color_scale[i] = (double *)malloc(NY*sizeof(double)); 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); 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++) 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) if (SAVE_TIME_SERIES)
{ {
wave_value = (long int)(phi[sample_left[0]][sample_left[1]]*1.0e16); wave_value = (long int)(phi[sample_left[0]][sample_left[1]]*1.0e16);
@ -760,8 +764,7 @@ void animation()
{ {
free(phi[i]); free(phi[i]);
free(psi[i]); free(psi[i]);
free(phi_tmp[i]); free(tmp[i]);
free(psi_tmp[i]);
free(total_energy[i]); free(total_energy[i]);
free(xy_in[i]); free(xy_in[i]);
free(color_scale[i]); free(color_scale[i]);

29
wave_comparison.c
View File

@ -232,7 +232,7 @@ double courant2, courantb2; /* Courant parameters squared */
/* animation part */ /* 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]) short int *xy_in[NX])
/* time step of field evolution */ /* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */ /* 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 */ /* evolve phi */
phi_out[i][j] = -y + 2*x + tcc[i][j]*delta - KAPPA*x - tgamma[i][j]*(x-y); 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++){ 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 */ /* evolve phi */
phi_out[i][j] = -y + 2*x + tcc[i][j]*delta - KAPPA*x - tgamma[i][j]*(x-y); 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; 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; 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; 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; 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; 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; 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 (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 */ /* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */ /* 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, psi, tmp, xy_in);
evolve_wave_half(phi_tmp, psi_tmp, phi, psi, 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() void animation()
{ {
double time, scale, energies[6], top_energy, bottom_energy; 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]; short int *xy_in[NX];
int i, j, s, counter = 0; int i, j, s, counter = 0;
@ -672,8 +663,7 @@ void animation()
{ {
phi[i] = (double *)malloc(NY*sizeof(double)); phi[i] = (double *)malloc(NY*sizeof(double));
psi[i] = (double *)malloc(NY*sizeof(double)); psi[i] = (double *)malloc(NY*sizeof(double));
phi_tmp[i] = (double *)malloc(NY*sizeof(double)); tmp[i] = (double *)malloc(NY*sizeof(double));
psi_tmp[i] = (double *)malloc(NY*sizeof(double));
xy_in[i] = (short int *)malloc(NY*sizeof(short int)); xy_in[i] = (short int *)malloc(NY*sizeof(short int));
} }
@ -763,7 +753,7 @@ void animation()
for (j=0; j<NVID; j++) 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); // 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(phi[i]);
free(psi[i]); free(psi[i]);
free(phi_tmp[i]); free(tmp[i]);
free(psi_tmp[i]);
free(xy_in[i]); free(xy_in[i]);
} }

26
wave_energy.c
View File

@ -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]) short int *xy_in[NX])
/* time step of field evolution */ /* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */ /* 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 */ /* evolve phi */
phi_out[i][j] = -y + 2*x + tcc[i][j]*delta - KAPPA*x - tgamma[i][j]*(x-y); 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; 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; 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; 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; 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 (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 */ /* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */ /* 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, psi, phi_tmp, psi_tmp, xy_in);
// evolve_wave_half_old(phi_tmp, psi_tmp, phi, psi, 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, psi, tmp, xy_in);
evolve_wave_half(phi_tmp, psi_tmp, phi, psi, 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() void animation()
{ {
double time, scale, energies[6], top_energy, bottom_energy; 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]; short int *xy_in[NX];
int i, j, s; int i, j, s;
@ -788,8 +782,7 @@ void animation()
{ {
phi[i] = (double *)malloc(NY*sizeof(double)); phi[i] = (double *)malloc(NY*sizeof(double));
psi[i] = (double *)malloc(NY*sizeof(double)); psi[i] = (double *)malloc(NY*sizeof(double));
phi_tmp[i] = (double *)malloc(NY*sizeof(double)); tmp[i] = (double *)malloc(NY*sizeof(double));
psi_tmp[i] = (double *)malloc(NY*sizeof(double));
xy_in[i] = (short int *)malloc(NY*sizeof(short int)); xy_in[i] = (short int *)malloc(NY*sizeof(short int));
} }
@ -852,7 +845,7 @@ void animation()
for (j=0; j<NVID; j++) 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); // 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(phi[i]);
free(psi[i]); free(psi[i]);
free(phi_tmp[i]); free(tmp[i]);
free(psi_tmp[i]);
free(xy_in[i]); free(xy_in[i]);
} }