restart;

# Halbwertszeiten
lambda_230 := ln(2)/7.538e4;
lambda_234 := ln(2)/2.455e5;
lambda_231 := ln(2)/3.276e4;

# Verhältnisse
ThU := Th_230/U_238 = 0.868; #+0.001; 
UU := U_234/U_238 = 1.0871; #+0.002;
UU_0 := U_234_0/U_238_0 = 1.14;
PaU := Pa_231/U_235 = 0.942; #+0.001;
PaU_0 := Pa_231_0/U_235_0 = 0;

Th_230 := solve(ThU,Th_230);
U_234 := solve(UU,U_234);
U_234_0 := solve(UU_0,U_234_0);
Pa_231 := solve(PaU,Pa_231);
Pa_231_0 := solve(PaU_0,Pa_231_0);

# 238U-234U
eq_UU := U_234/U_238 = U_234_0/U_238_0 * exp(-lambda_234*t);
t_UU := fsolve(eq_UU,t);

# TODO: check that
# 238U-230Th
eq_UTh := Th_230/U_234 = U_238/U_234 * (1-exp(-lambda_230*t)) + (1-U_238/U_234) * (lambda_230/(lambda_230-lambda_234)) * (1-exp(-(lambda_230-lambda_234)*t));
t_UTh := fsolve(eq_UTh,t);

# 231Pa-235U
eq_PaU := Pa_231/U_235 = Pa_231_0/U_235_0 * exp(-lambda_231*t)+1-exp(-lambda_231*t);
t_PaU := fsolve(eq_PaU,t);

# Alles auf einen Blick
# 238U-234U, 238U-230Th, 231Pa-235U
t_UU; t_UTh; t_PaU;