public static double Theta(double x, double a, double buf_max = 80, int k_max = 500) {
double sum, buf, xx, buf2;
sum += 2 * Math.Exp(-buf) * Math.Cos(2 * k * x);
} while ((buf < buf_max) && (k < k_max));
k = (int)Math.Round(x / (Math.PI));
xx = x - 2 * k * Math.PI;
sum = Math.Exp(-(a * xx * xx));
buf = (a * Math.Pow(xx + Math.PI * k, 2));
buf2 = (a * Math.Pow(xx - Math.PI * k, 2));
sum += Math.Exp(-buf) + Math.Exp(-buf2);
} while (((buf < buf_max) || (buf2 < buf_max)) && (k < k_max));
return Math.Sqrt(Math.PI * a) * sum;
private static double mz_ini(double a, int maxBuff = 80) {
mul = Math.Pow(Math.PI, 2) * a;
double sum = Math.Exp(-mul / 4);
buff = mul * Math.Pow(-2 * r + 0.5, 2);
sum += (-4 * r + 1) * Math.Exp(-buff);
buff = mul * Math.Pow(2 * r + 0.5, 2);
sum += (4 * r + 1) * Math.Exp(-buff);
} while (buff < maxBuff);
return sum * Math.PI * a * Math.Sqrt(Math.PI * a);
public static double[] CalcCoeffs(int n_cnt, int n_DCT, double a, int maxBuff = 80)
double[] work = new double[maxk_cnt];
double[] c = new double[maxk_cnt];
for (int m = 0; m <= n_DCT; m++)
work[m] = 1 / Math.Sqrt(Theta(m * h / 2, a));
for (int k = 0; k <= n_cnt; k++)
sum = (work[0] - work[n_DCT]) / 2;
sum = (work[0] + work[n_DCT]) / 2;
for (int m = 1; m < n_DCT; m++)
sum += work[m] * Math.Cos(k * h * m);
public static void Main()
double[] c_k = Utilities.CalcCoeffs(4, 499, 1, 80);