/* Program to write source file for LMODEL */
/* Based on program by                     */
/* John Vidale and Art Frankel, 1984       */
/* rewritten by Richard Stead, 7/7/87      */
/*      STYLE 0 - EXPLOSION                */
/*      STYLE 1 - STRIKE-SLIP EQ.          */
/*      STYLE 2 - DIP-SLIP EQ.             */



#include       <stdio.h>
#include       <math.h>
#include       <sys/file.h>

#define F       (float)
#define PI      3.1415926

void con();

main(ac,av) 
int ac; 
char **av;
{
  float dt;              /* time sampling rate */
  int   nt;              /* # of time pts in FD run */
  float h;               /* grid spacing, same units as dt, and srcvel */
  int   source;          /* inner radius of source box */
  int   style;           /* source mechanism */
  float trap1= F 0.0;    /* trapezoid rise */
  float trap2= F 0.0;    /* trapezoid sustain */
  float trap3= F 0.0;    /* trapezoid fall */
  float alpha= F 0.0;    /* gaussian half-width */
  float forcefac= F 0.0; /* line-force parameter for explosion adjustment */
  float srcden;          /* source region density */
  float srcvel;          /* source region velocity */
  float *y, *ts;
  float *xx;
  int sf, aint, space, source0, cnt, i, k;
  float sum, a1, r, conv;
  float x, z;
  float *pts;
  float *pxx, amp;
  int t1, t2, t3, limit;
  int atime;
  char sourcefile[64];
  float *srcout, *po;

  setpar(ac,av);
  mstpar("dt",      "f",&dt        );
  mstpar("nt",      "d",&nt        );
  mstpar("h",       "f",&h         );
  mstpar("source",  "d",&source    );
  mstpar("style",   "d",&style     );
  getpar("trap1",   "f",&trap1     );
  getpar("trap2",   "f",&trap2     );
  getpar("trap3",   "f",&trap3     );
  getpar("alpha",   "f",&alpha     );
  getpar("forcefac","f",&forcefac  );
  mstpar("source1", "s", sourcefile);
  mstpar("srcvel",  "f",&srcvel    );
  mstpar("srcden",  "f",&srcden    );
  endpar();

  sf = open(sourcefile,O_WRONLY|O_TRUNC|O_CREAT,0644);
  ts = (float *) malloc(4 * nt);

  if (trap1 == F 0.0 && trap2 == F 0.0 && trap3 == F 0.0) {
    /* gaussian time function */
    /* alpha<0, normalized one-sided gaussian */
    /* alpha>0, two-sided gaussian */

    aint = abs((int)(alpha * F 3.0));
    limit = 2 * aint + 1;
    xx = (float *) malloc(4 * (limit + 1));
    if (alpha >= F 0.0)
      for (i = -aint; i <= aint; i++)
        xx[i+aint] = F i * exp( F (-i * i) / (alpha * alpha));
    else {
      sum = F 0.0;
      for (i = -aint; i <= aint; i++) {
        xx[i+aint] = exp( F (-i * i) / (alpha * alpha));
        sum += xx[i+aint];
      }  
      for (i = 0; i <= limit; i++) xx[i] /= sum;
    }
  }  
  else {
    /* trapezoid time function */

    t1 = (int) (trap1 / dt);
    t2 = (int) (trap2 / dt);
    t3 = (int) (trap3 / dt);
    t2 += t1;
    limit = t2 + t3;
    xx = (float *) malloc(4 * (limit + 1));
    amp = (trap2 + F 0.5 * (trap1 + trap3)) / dt;
    pxx = xx;
    for (i = 0 ; i < t1; i++) *pxx++ = F i / ( F t1 * amp);
    for (      ; i < t2; i++) *pxx++ = F 1.0 / amp;
    for (i = t3; i > 0 ; i--) *pxx++ = F i / ( F t3 * amp);
  }
  y = (float *) malloc(4 * (limit + 2 + nt));

  /* source0 is innermost source ring radius */
  /* space is size of entire source file */
  source0 = source;
  space = 4 * source0 + 10;
  po = srcout = (float *) malloc(4 * nt * space);

  /* coefficients independent of r, theta, and t */
  if (style == 0)
    conv = sqrt( F 2.0 / srcvel) / PI;
  else
    conv = sqrt( F 2.0 / srcvel) / ( F 4.0 * PI * PI * srcden*srcvel*srcvel);

  /* loop over 4 rings, starting with outermost ring */
  cnt = 4;
  while (cnt--) {
    /* for each distance on ring */
    for (k=0; k<source+1; k++) {
      z = F source;
      x = F (source-k);
      r = sqrt(x * x + z * z);
      r *= h;
      /* FINDS TIME SERIES FOR EACH DISTANCE IN GIVEN RING */
      a1 = r / (srcvel * dt);
      atime = (int) a1;
      pts = ts;
      if (style == 0)
      {
        for (i = 0; i <= atime; i++) *pts++ = F 0.0;
        *pts++ = ( F 2.0 * conv / dt) * sqrt(( F 0.5 + F i - a1) / ( F 2.0*a1));
        i++;
        a1 *= a1;
        for ( ; i < nt; i++)
          *pts++ = conv / (dt * sqrt( F (i * i) - a1));
  
        /* Bad line??
         ts[i] = (2.*conv/dt)*sqrt(1/(i+a1))*(sqrt(i-a1+.5)-sqrt(i-a1-.5));
        */
      }
      else
      {
        for (i = 0; i <= atime; i++) *pts++ = F 0.0;
        *pts++ = ( F 2.0 * conv * F i / r)*sqrt(( F 0.5 + F i -a1)/( F 2.0*a1));
        i++;
        a1 *= a1;
        for ( ; i < nt; i++)
          *pts++ = conv / (r * sqrt( F 1.0 - a1 / F (i * i)));
      }

      /* CONVOLVE SOURCE FUNCTION */
      con(ts,xx,nt,limit,y);

      /* TRANSPOSE TO ORDER NEEDED IN PROGRAM */
      for (i = 0; i < nt; i++) po[i*space] = ts[i];
      po++;
    }
    source++;
  }
  write(sf,srcout,4*nt*space);
}

void
con(ts,x,nt,limit,y)
float *ts, *x;
register int nt;
int limit;
float *y;
{
  int i;
  register float *py, *pt, xi;
  register int j;

  py = y;
  for (i = 0; i < limit + nt + 1; i++) *py++ = F 0.0;
  for (i = 0; i < limit; i++) {
    py = y + i;
    xi = x[i];
    pt = ts;
    for (j = 0; j < nt; j++)
      *py++ += *pt++ * xi;
  }
  py = y;
  pt = ts;
  for (i = 0; i < nt; i++) *pt++ = *py++;
}