How can I obtain accurate phase information using fft from my multisine signal?

I have a sum-of-sine signal with different phases, I want to design a code which can reconstruct the original signal from the fft. I can do the power reconstruction. But my phase values in the reconstructed multisine is not matching my original signal. What am I doing wrong?

Here is my code so far:

    % code
inp = (repmat(3*sin(2*pi*1*(0:1/fs:1-1/fs)) + 1.5*sin(2*pi*31*(0:1/fs:1-1/fs) + pi/2),200,1))';
ph = (fft(inp));
ph(abs(ph)<1e-6) = 0;
ang = angle(ph);
inpspk = fft(inp,fs)./fs;
inpspk = inpspk.*conj(inpspk);
pow = inpspk;
inpspk = mean(inpspk,2);
mysig = zeros(1000,200);
for j = 1:Epochs
    for i = 2:498
        mysig(:,j) = mysig(:,j) + (sqrt(pow(i+1,j)) * sin(2 * pi * i * (0:1/fs:1-1/fs) + ang(i+1,j)))';
    end
end
ppt = zeros(1,100);
for j = 1:length(ppt)
    ppt(j) = inpspk(j + 1)/sum(inpspk(2:101)) * 100;
end
%Plot input spectrum
figure
stem (1:length(ppt),ppt, 'LineWidth', 2);
xlabel('Frequency (Hz)','Fontsize',24)
set(gca,'FontSize',24);
grid on;
recon = mysig;
figure
plot(recon(:,1));
ph = (fft(recon));
ph(abs(ph)<1e-6) = 0;
ang2 = angle(ph);
recon = fft(recon,fs) ./ fs;
recon = recon.*conj(recon);
recon = mean(recon,2);
ppt = zeros(1,100);
for j = 1:length(ppt)
    ppt(j) = recon(j + 1)/sum(recon(2:101)) * 100;
end
%Plot output spectrum
figure
stem (1:length(ppt),ppt, 'LineWidth', 2);
xlabel('Frequency (Hz)','Fontsize',24)
set(gca,'FontSize',24);
grid on;
figure
plot(ang(:,1)/pi);
hold on
plot(ang2(:,1)/pi);

[KALYAN ACHARJYA]

Format your code. Read here

[Nirvik Sinha]

Thank you for the advice. I have done it. Please help me with the problem.