How to perform sine generation (in C)? I am interested in how should a listing for a program look like that will generate a sine of f=50hz, using the next arguments for any moment in time.

float phis=162*PI/180;//initial phase of signal
//////////////////////////////////////////////////////////////////////// 
//  generates a sinusoid of dwLength in dwData 
//      and maintains phis value updated
////////////////////////////////////////////////////////////////////////
void rc(short* dwDataOut, DWORD dwLength){
float omegaf=50/(float)dwLength/PI*90/100;
DWORD dwi;
//use phis
    for (dwi=0;dwi&lt;dwLength;dwi++)
        dwDataOut[dwi]=4450*sin(dwi*omegaf+phis)+200;

//update wave phase
    phis-=PI*(1-50/(float)dwLength)+PI/2;
    while(phis&lt;0)
        phis+=2*PI;
}


The program I came up works correctly with a small marj of error only for f=50hz, because right from the start it appears with a small error in getting phis coefficient.

By applying this procedure in a program I'll make use of a cleaner signal from a sound wave, and obtain a signal that will be overlapped as antiphase signal to a more complex 50hz wave, for longer period of time, and so eliminating this harmonic. I’m interested in a more mathematical method to find the value for phis, and probably the arguments of sine function.

(is this function resulting in a true output for f=50hz? what are the errors that can be adjusted?)
I think your recalculation of phis in the function is not correct.

If I do this in scilab:

//25631
N = 100;
phis1=162*%pi/180;//initial phase of signal

[data1,phis2] = rc(N,phis1);
[data2,phis3] = rc(N,phis2);
[data3,phis4] = rc(N,phis3);

clf
plot([data1 data2 data3])


where rc is 

//////////////////////////////////////////////////////////////////////// 
//  generates a sinusoid of dwLength in dwData 
//      and maintains phis value updated
////////////////////////////////////////////////////////////////////////
function [dwDataOut,phis] = rc(dwLength, phis)
omegaf=50/dwLength/%pi*90/100;

dwDataOut = zeros(1,dwLength);

for dwi=1:dwLength,
    dwDataOut(dwi) = 4450*sin(dwi*omegaf+phis)+200;
end

//update wave phase
//WAS: phis = phis - %pi*(1-50/dwLength)+%pi/2;
phis = dwLength*omegaf+phis;
while(phis&lt;0)
    phis = phis + 2*%pi;
end
endfunction


Note specifically the line phis = dwLength*omegaf+phis;

Then the plot is as below, which seems to be what you're after.