% spall_fd_vs_exact.m -- Kabuk kirilmasi: acik sonlu fark semasi vs KESIN cozum
% Dogrusal elastik dalgada yansiyan cekme gelen tepeyi ASAMAZ; kalin duvar +
% kisa darbede kesin sonuc sigma_t/sigma_0 = 1.
function spall_fd_vs_exact()
  E = 30e9; nu = 0.2; rho = 2400;
  t = 1.0;  Pr = 5e6; ir = 250;
  td = 2*ir/Pr;
  M  = E*(1-nu)/((1+nu)*(1-2*nu));
  c  = sqrt(M/rho); T = t/c;
  load = @(s) (s >= 0 & s <= td) .* Pr .* (1 - s/td);

  % --- 1) Acik sonlu fark (orijinal spallCheck.m semasi) ---
  nel = 80; he = t/nel; dt = 0.5*he/c;
  mN = rho*he*ones(nel+1,1); mN([1 end]) = mN([1 end])/2;
  u = zeros(nel+1,1); v = zeros(nel+1,1); smin_fd = 0;
  for tt = 0:dt:(td + 4*T)
    s = M*(u(1:end-1) - u(2:end))/he;
    f = zeros(nel+1,1); f(1) = load(tt);
    f(1:end-1) = f(1:end-1) - s;  f(2:end) = f(2:end) + s;
    v = v + dt*f./mN;  u = u + dt*v;
    smin_fd = min(smin_fd, min(s));
  end

  % --- 2) KESIN cozum (d'Alembert) ---
  x = linspace(0, t, 81); smin_ex = 0;
  for tau = linspace(0, td + 6*T, 3001)
    sig = F(tau - x/c, load, T) - F(tau + x/c - 2*T, load, T);
    smin_ex = min(smin_ex, min(sig));
  end
  fprintf('dalga hizi c = %.0f m/s   td = %.3f ms   T = t/c = %.1f us\n', c, td*1e3, T*1e6);
  fprintf('Sonlu fark : sigma_t/sigma_0 = %.3f\n', -smin_fd/Pr);
  fprintf('Kesin cozum: sigma_t/sigma_0 = %.3f   (kuram: 1.000)\n', -smin_ex/Pr);
  fprintf('Sonlu fark asimi = %%%.1f\n', (smin_fd/smin_ex - 1)*100);
end

function out = F(tau, load, T)
  out = zeros(size(tau));
  for kk = 0:(floor(max(tau)/(2*T)) + 1)
    out = out + load(tau - 2*kk*T);
  end
end