import pymatbridge
ip = get_ipython()
pymatbridge.load_ipython_extension(ip)

%%matlab
path(path,'/Users/rjl/SMM/all_M_files')

%%matlab

% p5.m - repetition of p4.m via FFT
%        For complex v, delete "real" commands.

% Differentiation of a hat function:
  N = 24; h = 2*pi/N; x = h*(1:N)';
  v = max(0,1-abs(x-pi)/2); v_hat = fft(v);
  w_hat = 1i*[0:N/2-1 0 -N/2+1:-1]' .* v_hat;
  w = real(ifft(w_hat)); clf
  subplot(3,2,1), plot(x,v,'.-','markersize',13)
  axis([0 2*pi -.5 1.5]), grid on, title('function')
  subplot(3,2,2), plot(x,w,'.-','markersize',13)
  axis([0 2*pi -1 1]), grid on, title('spectral derivative')

% Differentiation of exp(sin(x)):
  v = exp(sin(x)); vprime = cos(x).*v;
  v_hat = fft(v);
  w_hat = 1i*[0:N/2-1 0 -N/2+1:-1]' .* v_hat;
  w = real(ifft(w_hat));
  subplot(3,2,3), plot(x,v,'.-','markersize',13)
  axis([0 2*pi 0 3]), grid on
  subplot(3,2,4), plot(x,w,'.-','markersize',13)
  axis([0 2*pi -2 2]), grid on
  error = norm(w-vprime,inf);
  text(2.2,1.4,['max error = ' num2str(error)])


%%matlab
N = 16;
x = linspace(2*pi/N,2*pi,N);
ik = 1i*[0:N/2 -N/2+1:-1];   % i * wave number vector (matlab ordering)
ik2 = ik.*ik;                    % multiplication factor for second derivative

u = exp(cos(x));
u_hat = fft(u);
v_hat = ik2 .* u_hat;
v = real(ifft(v_hat));     % imaginary part should be at machine precision level

error = v - (sin(x).^2 - cos(x)) .* exp(cos(x));
norm(error,inf)

%%matlab

N = 16;
x = linspace(2*pi/N,2*pi,N);
f = (sin(x).^2 - cos(x)) .* exp(cos(x));
f_hat = fft(f);

ik = 1i*[0:N/2 -N/2+1:-1];   % i * wave number vector (matlab ordering)
ik2 = ik.*ik;         % multiplication factor for second derivative
ii = find(ik ~= 0);   % indices where ik is nonzero
ik2inverse = ik2;     % initialize zeros in same locations as in ik2
ik2inverse(ii) = 1./ik2(ii);   % multiplier factor to solve u'' = f

u_hat = ik2inverse .* f_hat;
u = real(ifft(u_hat));   % imaginary parts should be roundoff level

%%matlab
plot(x,u,'b-o')
hold on
v = exp(cos(x));
plot(x,v,'r-o')

%%matlab
u2 = u + v(1)-u(1);
norm(u2 - v, inf)

from scipy import fft,ifft
N = 16;
x = linspace(2*pi/N,2*pi,N)

ik = 1j*hstack((range(0,N/2+1), range(-N/2+1,0)));   # i * wave number vector (matlab ordering)
ik2 = ik*ik;          # multiplication factor for second derivative

u = exp(cos(x))
u_hat = fft(u)
v_hat = ik2 * u_hat
v = real(ifft(v_hat))     # imaginary part should be at machine precision level

error = v - (sin(x)**2 - cos(x)) * exp(cos(x))
norm(error,inf)




N = 16;
x = linspace(2*pi/N,2*pi,N)
f = (sin(x)**2 - cos(x)) * exp(cos(x))
f_hat = fft(f)
ik = 1j*hstack((range(0,N/2+1), range(-N/2+1,0)));   # i * wave number vector (matlab ordering)
ik2 = ik*ik;          # multiplication factor for second derivative
ik2inverse = where(ik2 != 0, 1./ik2, 0.)
u_hat = ik2inverse * f_hat;
u = real(ifft(u_hat))

plot(x,u,'b-o')
v = exp(cos(x));
plot(x,v,'r-o')



u2 = u + v[0]-u[0]
norm(u2 - v, inf)