// propagator.r -- This routine constructs a propagator at the
//                 specified time and plots the real part.  The
//                 propagator, the x axis vector, the propagator
//                 type, and the time are returned as a list.

// Argument list:  dispreln -- This is an integer with the range
//                      1-3.  1 implies omega = |k|.  2 implies
//                      omega = k*k/2.  3 implies omega = sqrt(1 + k*k).
//                 time -- This is the time at which the propagator
//                      is evaluated.  Values in the range 0-50
//                      make sense, since the x domain is [-50, 50].
//                      A wave moving at the speed of light (c = 1)
//                      will move from the origin to the edge of the
//                      domain at time = 50.

// Suggestion: Run first with time = 0 to verify Dirac delta function
//             behavior.  Then try all dispersion relations for time = 2.
//             Finally, gradually increase the time to see how the
//             propagator evolves for each dispersion relation.

// Theory: The propagator is defined as K(x,t) = integral from -inf to
//         inf of exp[i(k*x - omega*t)] dk/(2*pi).  This is an ill-
//         defined integral, but an approximation can be made by tapering
//         off the value of the integrand for large values of k.  This
//         is done here with the tapering function (1 + cos(pi*k/kmax))/2,
//         where the integral is numerically integrated over [-kmax,kmax]
//         rather than [-inf,inf].

// This is the main routine which calculates the propagator and plots it.
propagator = function(dispreln,time) {
	if (nargs != 2) {
		error("Usage: propagator(dispreln,time)");
	}
// construct some preliminary stuff
	nx = 250;
	dx = .2;
	x = [-nx:nx]'*dx;
	dk = 6.28318/(2*nx*dx);
	k = [-nx+1:nx]*dk;
	t = 0*x + time;
// compute the frequency vector
	if (dispreln < 1 || dispreln > 3) {
		error("propagator: allowed dispersion relations are 1-3");
	}
	if (dispreln == 1) {
		om = abs(k);
		label = "omega = |k|";
	}
	if (dispreln == 2) {
		om = k.*k/2;
		label = "omega = k*k/2";
	}
	if (dispreln == 3) {
		om = sqrt(1 + k.*k);
		label = "omega = sqrt(1 + k*k)";
	}
// compute the taper factor
	factor = 0*x + 6.28318/(2*nx*dk);
	taper = (1 + cos(factor*k))/2;
// compute the integrand
	integrand = exp(1i*(x*k - t*om)).*taper;
// integrate
	kernel = sum(integrand')'*dk/6.28318;
// do the plot
	pgopen();
	pgclear();
	pgwindow(1,0,0,15,10);
	vert = 1/(2*dx);
	horiz = nx*dx;
	pgaxes(1,-horiz,horiz,11,"x",-vert,vert,11,"real(kernel)");
	pgline(1,5,x,real(kernel));
	sprintf(tlabel,"Time: %4.1f; Dispersion relation: %s",time,label);
	pglabel(1,-.9*horiz,.8*vert,tlabel);
	pgpause();
	pgclose();
// return the stuff
	return <<propagator = kernel; x = x; dispreln = label; time = time>>;
}