Get spectrum from autocorrelation function using fft()

I would like to observe Fourier Spectrum, got through Autocorrelation function. This is my code:

clear;clc;

%% generate x-axe
first_step=0;
step_t=0.01;
last_step=1;
t=first_step:step_t:last_step;

%% signal and fourier transform
y=  2*sin(2*pi*60*t) ;
figure(1)
plot(t,y);title('signal')
fourier=fft(y );
N_=length(fourier);
f_ = (1 / step_t)  *   (   0:    (N_/2)   )    /    N_;
a=(fourier.*conj(fourier))/(N_*N_);  % Spectral Density get half
a = a(1:N_ /2+1);
a(2:end-1) = 4*a(2:end-1);
figure(2);
plot(f_,a);title('Spectrum Power through signal');

%% through AutoCorrelation function
tau=0.01;
y1=  2*sin(2*pi*60*(t+tau)) ;
Y=y1.*y;
AutoCorr_func = cumtrapz(Y,t);

figure(3);
plot(t,AutoCorr_func); title('AutoCorrelation function');
fourier=fft(AutoCorr_func);
N_=length(fourier);
f_ = (1 / step_t)  *   (   0:    (N_/2)   )    /    N_;
a=fourier.*conj(fourier) ; % Spectral Density
a = a(1:N_ /2+1);          % spectrum is even, so get half
a(2:end-1) =  a(2:end-1);  % multiply by 4, since spectrum is square of amplitude, and we have two even halves

figure(4)
plot(f_ ,a );
title('Spectrum through AutoCorrelation funciton');

so a few questions:

1) why I get so strang ACF? It has to be periodic, since y and y1 are strongly periodic. Is sonething wrong with it? maybe here

Y=y1.*y;
AutoCorr_func = cumtrapz(Y,t);

but Acf is an integral and I am trying to get an array numerically.

2) My tries to get power density through Acf are the same as in the simple case, without using Acf, but Figure 2 gives correct resultL frequency = 42 and power is square of amplitude, in my case it should be 4, but I have around 2. What is the problem...?

3) How to get correct power density through Acf, what am I doing wrong?

Any help would be greatly appreciate! Thank you in advance!