% Interfase metal-superconductor
% Adimensionalizo:
% - longitud en unidades de lambda_L
% - corriente en unidades de j_metal

d = 5;

X = linspace(-d, d, 201);
Y = linspace(-0.5*d, 1.5*d, 401);

[XX, YY] = meshgrid (X, Y);

lambda = @(n) 1 ./ sqrt(1 + (n * pi / d).^2);

c_x = @(n) 2 * (-1).^(n + 1) .* lambda(n) ./ n * d / pi;
c_y = @(n) 2 * (-1).^(n + 1) .* lambda(n).^2;

% n = 0
jx = zeros(size(XX));
jy = d/sinh(d) *  cosh(XX);

n_max = 100;

for (n = 1:n_max)
    jx += c_x(n) * exp(-YY / lambda(n)) .* sin(n * pi * XX / d);
    jy += c_y(n) * exp(-YY / lambda(n)) .* cos(n * pi * XX / d);

endfor

jx = ifelse(YY > 0, jx, zeros(size(XX)));
jy = ifelse(YY > 0, jy, ones(size(XX)));

figure(1);
pcolor(XX,YY,jx)
shading interp
axis("equal")
colorbar()

figure(2);
pcolor(XX,YY,jy)
shading interp
axis("equal")
colorbar()


f = @(x,y) interp2(XX,YY,jx, x, y)./interp2(XX,YY,jy, x, y);

Xi = linspace(-0.99*d, 0.99*d, 15);
Yi = linspace(-0.05*d, 1.45*d, 101);

Xl = lsode(f, Xi, Yi);

figure(3); clf()
colormap('hot');

pcolor(YY',XX', abs(jx + i*jy)')
shading interp
axis("equal")
xlabel("y / lambda_L")
ylabel("x / lambda_L")

cb = colorbar()
cb_ylab = get(cb,"ylabel")
set(cb_ylab,"string","|j_S| / j_M ")

hold on
plot(Yi,Xl,'-g')
axis("equal")