#include <stdio.h>
#include <math.h>
/*-----------------------------------------------------
    IIR digital filter to convert acceleration to
    displacement record of a mechanical seismograph
    Y(n) = gn * (X(n) + h[0] * X(n-1) + h[1] * X(n-2)
             - h[2] * Y(n-1) - h[3] * Y(n-2)
    return code =0 : no error
                 =1 : t0 <= 0
                 =2 : dt <= 0
  -----------------------------------------------------*/
int filt_acc_to_seismogram(
       double t0,         /* period of seismograph */
       double h_damp,     /* damping constant */
       double dt,         /* SAMPLING INTERVAL */ 
       int    *m_filt,    /* DEGREE OF FILTER */
       double *gn_filt,   /* GAIN FACTOR */
       double h_filt[])   /* FILTER CONSTANTS */
{
   double  pi= M_PI;
   double  c,
           w0;

   if(t0 <= 0.0) return(1);
   if(dt <= 0.0) return(2);

   c=(2.0 * pi / t0)/tan(pi * dt / t0);
   w0=2.0 * pi /(c * t0);
   *gn_filt =1.0/(1.0+2.0 * h_damp * w0 + w0 *w0) 
             * (dt/2.0) * (dt/2.0);
   h_filt[0]= 2.0;
   h_filt[1]= 1.0;
   h_filt[2]=(-2.0+2.0 * w0 * w0)
            /(1.0+2.0 * h_damp * w0 + w0 * w0);
   h_filt[3]=(1.0-2.0 * h_damp * w0 + w0 * w0)/
             (1.0+2.0 * h_damp * w0 + w0 * w0);
   *m_filt =1;
   return(0);
}