单质点阻尼系的地震反应(线性加速度法-C#代码)

        /// <summary>
        /// 求解固有频率为f的单自由度在给定加速度时程ath激励下的加速度时程响应
        /// </summary>
        /// <param name="h">阻尼比</param>
        /// <param name="f">单自由度系统固有频率</param>
        /// <param name="dt">加速度时程采样时间间隔</param>
        /// <param name="ath">加速度时程</param>
        /// <returns></returns>
        public static void  SdofResponse(double h, double f, double dt, double[] ath,ref double[] acc,ref double[] vel,ref double[] dis)
        {
            int nSample = ath.Length;
             //acc = new double[nSample];
             //vel = new double[nSample];
             //dis = new double[nSample];
            double w = 2 * Math.PI * f;
            double w2 = w * w;
            double hw = h * w;
            double wd = w * Math.Sqrt(1 - h * h);
            double wdt = wd * dt;
            double E = Math.Exp(-hw * dt);
            double cwdt = Math.Cos(wdt);
            double swdt = Math.Sin(wdt);
            double a11 = E * (cwdt + hw * swdt / wd);
            double a12 = E * swdt / wd;
            double a21 = -E * w2 * swdt / wd;
            double a22 = E * (cwdt - hw * swdt / wd);

            double s1 = 2 * h / w2 / w / dt;
            double s2 = (1 - 2 * h * h) / w2 / wdt;
            double s3 = 1 / w2 / dt;
            double s4 = h / w / wdt;
            double b11 = E * ((1 / w2 + s1) * cwdt + (h / w / wd - s2) * swdt) - s1;
            double b12 = E * (-s1 * cwdt + s2 * swdt) - 1 / w2 + s1;
            double b21 = E * (-s3 * cwdt - (s4 + 1 / wd) * swdt) + s3;
            double b22 = E * (s3 * cwdt + s4 * swdt) - s3;

            acc[0] = 2 * hw * ath[0] * dt;
            vel[0] = -ath[0] * dt;
            dis[0] = 0;
            //acc[0] = 0;
            //vel[0] = 0;
            //dis[0] = 0;

            for (int i = 1; i < nSample; i++)
            {
                dis[i] = a11 * dis[i - 1] + a12 * vel[i - 1] + b11 * ath[i - 1] + b12 * ath[i];
                vel[i] = a21 * dis[i - 1] + a22 * vel[i - 1] + b21 * ath[i - 1] + b22 * ath[i];
                acc[i] = -2 * hw * vel[i] - w2 * dis[i];
            }
        }