Yes, there is a correct way and it's simple, too.
By definition, the autocorrelation of a stationary process $X_t$ at lag $dt$ is the correlation between $X_t$ and $X_{t+dt}$. Suppose you have observations of this process $x_{t_0}, x_{t_0+dt}, x_{t_0+2dt}, \ldots, x_{t_0+k_0dt}$ at lag $dt$, another set of observations in a non-overlapping time interval $x_{t_1}, x_{t_1+dt}, x_{t_1+2dt}, \ldots, x_{t_1+k_1dt}$ at lag $dt$ for $t_1 \gt t_0+d_0dt$, and in general you have contiguous observations of samples $x_{t_i}, x_{t_i+dt}, x_{t_i+2dt}, \ldots, x_{t_i+k_idt}$, $i=0, 1, \ldots$ for non-overlapping time intervals. Then the correlation coefficient of the ordered pairs
$$\{(x_{t_i+jdt}, x_{t_i+(j+1)dt})\}$$
for $i=0, 1, \ldots$ and $j=0, 1, k_i-1$ estimates the autocorrelation of $x_t$ at lag $dt$. Compute the standard errors of the correlation exactly as you would compute the standard error for the correlation of any bivariate data set $\{(x_k, y_k)\}$.
The difference between this approach and the one proposed in the question is that pairs spanning two sequences, $(x_{t_j+k_jdt}, x_{t_{j+1}})$, are not included in the calculation. Intuitively they should not be, because in general the time interval between these pairs is not equal to $dt$ and therefore such pairs do not provide direct information about the correlation at lag $dt$.
If you make your harmonic signal to have a zero mean you'd see the auto correlation have more negative values.
Since your signal is almost always positive and the correlation is a sum of the values it makes sense it will be almost always positive.
Best Answer
Have a look here. General idea is to generate normal random variables with give covariance matrix.
If Z is a vector of uncorrelated Gaussian random variables, and C is a square root of the given covariance matrix, then target variable will be mean+CZ.
To obtain C You may use Cholesky decomposition.
If You may obtain Gaussian variable with a given covariance matrix, then You may convert it to differently distributed variable using inverse CDF method.
Important to notice, that CDF of any target vector is uniformly distributed, than You just map your CDF to a desired. Finally You have random variable from any distribution wich has given covariance matrix.