/* 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 - P-WAVE EXPLOSION	   */
/*	STYLE 1 - S-WAVE EXPLOSION	   */
/*	STYLE 2 - STRIKE-SLIP EQ.	   */
/*	STYLE 3 - DIP-SLIP EQ.		   */
/*	STYLE 4 - 45 DIP-SLIP EQ.	   */


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

#define F       (float)
#define PI	(float) 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 srcvels */
  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 srcpvel;         /* source region p-velocity */
  float srcsvel;         /* source region s-velocity */
  char line[256];
  float *y, *tssu, *tssw, *tspu, *tspw, *tsu, *tsw;
  float *xx;
  int sf, aint, space, source0, cnt, i, j, k;
  float sum, a1, r, cons, conp, s2, c2, fi;
  float x, z, *ps, *px, *pz, sx, sz, ss;
  float holdu;
  float *ppu, *ppw, *psu, *psw;
  float *pxx, amp;
  int t1, t2, t3, limit;
  int atime, btime;
  float f1, f2, sT, cT;
  char sourcefile[64];
  char tempfile[64];
  float loopc1, loopc2, loopc3, loopc4, loopc5, loopc6, loopc7, loopc8;

  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("tempfile","s", tempfile  );
  mstpar("srcpvel", "f",&srcpvel   );
  mstpar("srcsvel", "f",&srcsvel   );
  mstpar("srcden",  "f",&srcden    );
  endpar();

  sf = open(tempfile,O_WRONLY|O_TRUNC|O_CREAT,0644);
        tssu = (float *) malloc(4 * nt);
        tssw = (float *) malloc(4 * nt);
  tsu = tspu = (float *) malloc(4 * nt);
  tsw = tspw = (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 = 32 * source0 + 48;

  /* coefficients independent of r, theta, and t */
  if (style == 0)
    conp = sqrt( F 2.0 / srcpvel) / PI;
  else
    conp = sqrt( F 2.0 / srcpvel) / ( F 4.0 * PI * PI * srcden);
  if (srcsvel != F 0.0)
    cons = sqrt( F 2.0 / srcsvel) / ( F 4.0 * PI * PI * srcden);
  else
    cons = F 0.0;

  /* loop over 4 rings, starting with outermost ring */
  cnt = 4;
  while (cnt--) {
    /* LOOPS OVER 4 SIDES OF EACH RING */
    px = &x; pz = &z; sx = sz = F 1.0;
    for (j = 0; j < 4; j++) {
      for (k = 0; k < 2 * source; k++) {
        *px = F (k - source) * sx;
        *pz = F source * sz;
        r = sqrt(x * x + z * z);
        sT = x / r;
        cT = z / r;
        r *= h;
        s2 = sT * sT;
        c2 = cT * cT;
        atime = (int) (r / (srcpvel * dt));
        btime = (int) (r / (srcsvel * dt));
        ppu = tspu;
        ppw = tspw;
        psu = tssu;
        psw = tssw;

        /* FINDS TIME SERIES FOR EACH GRID POINT ON SIDE */
        switch (style) {

        case 0:			/* EXPLOSION */
          for (i = 2; i <= atime; i++) {
            *ppu++ = F 0.0;
            *ppw++ = F 0.0;
          }

          /* RIG FIRST POINT AT SINGULARITY */
          a1 = conp / (dt * srcpvel);
          f1 = ( F i + F 0.5) * srcpvel * dt / r;
          f2 = sqrt(f1 * f1 - F 1.0);
          holdu = F 0.5 * (f1 * f2 + log(f1 + f2)) - c2 * f1 * f2;
          *ppu++ = (( F 1.0 - forcefac) * sT * f2 + forcefac * holdu) * a1;
          *ppw++ = ((forcefac - F 1.0) - forcefac * sT * f1) * cT * f2 * a1;

          /* old first point:
          tspw[i-2]=a1*sqrt(srcpvel*(dt+2.*(i*dt-r/srcpvel))/r);
          tspu[i-2]= sT*tspw[i-2];
          tspw[i-2] *= -cT;
          holdw=a1*i*dt*srcpvel/r*sqrt(srcpvel*(dt+2.*(i*dt-r/srcpvel))/r);
          holdu= s2*holdw;
          holdw *= -sT*cT;
          tspu[i-2] = (1-forcefac)*tspu[i-2] + forcefac*holdu;
          tspw[i-2] = (1-forcefac)*tspw[i-2] + forcefac*holdw;
          */

          i++;

          loopc1 = r / (dt * srcpvel);
          loopc2 = loopc1 * loopc1;
          loopc3 = ( F 1.0 - forcefac) * sT;
          loopc4 = (forcefac - F 1.0) * cT;
          loopc5 = forcefac * (s2 - c2) / loopc1;
          loopc6 = -forcefac * cT * (sT + sT) / loopc1;
          loopc7 = forcefac * c2 * loopc1;
          loopc8 = forcefac * cT * sT * loopc1;
          for ( ; i < nt + 1; i++) {
            fi = F i;
            *ppw   = conp / (r * sqrt( F 1.0 - loopc2/ (fi*fi)));
            *ppu++ = *ppw * (loopc3 + loopc5 * fi + loopc7 / fi);
            *ppw++ *= loopc4 + loopc6 * fi + loopc8 / fi;
          }

          for (i = 2; i < nt + 1; i++) {
            *psu++ = F 0.0;
            *psw++ = F 0.0;
          }
        break;

        case 1:			/* S-WAVE "EXPLOSION" */
          for (i = 2; i < nt + 1; i++) {
            *ppu++ = F 0.0;
            *ppw++ = F 0.0;
          }

          for( i = 2; i <= btime; i++) {
            *psu++ = F 0.0;
            *psw++ = F 0.0;
          }

          f1 = ( F i + F 0.5) * srcsvel * dt / r;
          f2 = sqrt(f1 * f1 - F 1.0);
          *psw    = F i * cons * f2 / r;
          *psu++  = cT * *psw;
          *psw++ *= sT;
          i++;

          loopc1 = r / (dt * srcsvel);
          loopc1 *= loopc1;
          loopc2 = cons / r;
          for( ; i < nt + 1; i++) {
            *psw    = loopc2 / sqrt( F 1.0 - loopc1 / F (i*i));
            *psu++  = cT * *psw;
            *psw++ *= sT;
          }
        break;

        case 2:			/* STRIKE-SLIP EQ */
          for(i = 2; i <= atime; i++) {
            *ppu++ = F 0.0;
            *ppw++ = F 0.0;
          }

          a1 = F i * dt / r;
          a1 *= a1 * a1 * conp;
          f1 = ( F i + F 0.5) * srcpvel * dt / r;
          f2 = sqrt(f1 * f1 - F 1.0);
          *ppw    = a1 * f2 / dt;
          *ppu++  =  s2 * sT * *ppw;
          *ppw++ *= -s2 * cT;
          i++;

          loopc1 = r / (dt * srcpvel);
          loopc1 *= loopc1;
          loopc2 = conp * dt * dt / (r * r * r);
          loopc3 = sT * (s2 - F 3.0 * c2);
          loopc4 = cT * (c2 - F 3.0 * s2);
          loopc5 = loopc1 * sT *   F 3.0 * c2;
          loopc6 = loopc1 * cT * ( F 2.0 * s2 - c2);
          for ( ; i < nt + 1; i++) {
            fi = F (i * i);
            *ppw    = loopc2 * fi / sqrt( F 1.0 - loopc1 / fi);
            *ppu++  = *ppw * (loopc3 + loopc5 / fi);
            *ppw++ *=         loopc4 + loopc6 / fi;
          }

          for (i = 2; i <= btime; i++) {
            *psu++ = F 0.0;
            *psw++ = F 0.0;
          }

          a1 = F i * dt / r;
          a1 *= a1 * a1 * cons;
          f1 = ( F i + F 0.5) * srcsvel * dt / r;
          f2 = sqrt(f1 * f1 - F 1.0);
          *psw    = a1 * f2 / dt;
          *psu++  = c2 * sT * *psw;
          *psw++ *= s2 * cT;
          i++;

          loopc1 = r / (dt * srcsvel);
          loopc1 *= loopc1;
          loopc2 = cons * dt * dt / (r * r * r);
          loopc3 = sT * (3.0 * c2 - s2);
          loopc4 = cT * (3.0 * s2 - c2);
          loopc5 = loopc1 * sT * (s2 - 2.0 * c2);
          loopc6 = loopc1 * cT * (c2 - 2.0 * s2);
          for ( ; i < nt + 1; i++) {
            fi = F (i * i);
            *psw    = loopc2 * fi / sqrt( F 1.0 - loopc1 / fi);
            *psu++  = *psw * (loopc3 + loopc5 / fi);
            *psw++ *=         loopc4 + loopc6 / fi;
          }
          break;

       case 3:			/* DIP-SLIP EQ */
          for (i = 2; i <= atime; i++) {
            *ppu++ = F 0.0;
            *ppw++ = F 0.0;
          }

          a1 = F i * dt / r;
          a1 *= a1 * a1 * conp;
          f1 = (F i + F 0.5) * srcpvel * dt / r;
          f2 = sqrt(f1 * f1 - F 1.0);
          *ppw    = a1 * f2 / dt;
          *ppu++  = F  2.0 * s2 * cT * *ppw;
          *ppw++ *= F -2.0 * c2 * sT;
          i++;

          loopc1 = r / (dt * srcpvel);
          loopc1 *= loopc1;
          loopc2 = conp * dt * dt / (r * r * r);
          loopc3 = cT * ( F 6.0 * s2 - F 2.0 * c2);
          loopc4 = sT * ( F 2.0 * s2 - F 6.0 * c2);
          loopc5 = loopc1 * cT * ( F 2.0 * c2 - F 4.0 * s2);
          loopc6 = loopc1 * sT * ( F 4.0 * c2 - F 2.0 * s2);
          for ( ; i < nt + 1; i++) {
            fi = F (i * i); 
            *ppw    = loopc2 * fi / sqrt( F 1.0 - loopc1 / fi);
            *ppu++  = *ppw * (loopc3 + loopc5 / fi);
            *ppw++ *=         loopc4 + loopc6 / fi;
          }

          for (i = 2; i <= btime; i++) {
            *psu++ = F 0.0;
            *psw++ = F 0.0;
          }

          a1 = F i * dt / r;
          a1 *= a1 * a1 * cons;
          f1 = ( F i + F 0.5) * srcsvel * dt / r;
          f2 = sqrt(f1 * f1 - F 1.0);
          *psw    = a1 * f2 / dt;
          *psu++  = (c2 * cT - cT * s2) * *psw;
          *psw++ *= (c2 * sT - sT * s2);
          i++;

          loopc1 = r / (dt * srcsvel);
          loopc1 *= loopc1;
          loopc2 = cons * dt * dt / (r * r * r);
          loopc3 = cT * ( F 2.0 * c2 - F 6.0 * s2);
          loopc4 = sT * ( F 6.0 * c2 - F 2.0 * s2);
          loopc5 = loopc1 * cT * ( F  5.0 * s2 - c2);
          loopc6 = loopc1 * sT * ( F -5.0 * c2 + s2);
          for ( ; i < nt + 1; i++) {
            fi = F (i * i);
            *psw    = loopc2 * fi / sqrt( F 1.0 - loopc1 / fi);
            *psu++  = *psw * (loopc3 + loopc5 / fi);
            *psw++ *=         loopc4 + loopc6 / fi;
          }
        break;

        case 4:			/* 45 DEGREE DIP-SLIP EQ */
          for (i = 2; i <= atime; i++) {
            *ppu++ = F 0.0;
            *ppw++ = F 0.0;
          }

          a1 = F i * dt / r;
          a1 *= a1 * a1 * conp;
          f1 = ( F i + F 0.5) * srcpvel * dt / r;
          f2 = sqrt(f1 * f1 - F 1.0);
          *ppw    = a1 * f2 / dt;
          *ppu++  = (-s2 * sT + F 2.0 * sT * c2) * *ppw;
          *ppw++ *= ( s2 * cT - F 2.0 * cT * c2);
          i++;

          loopc1 = r / (dt * srcpvel);
          loopc1 *= loopc1;
          loopc2 = conp * dt * dt / (r * r * r);
          loopc3 = sT * ( F 9.0 * c2 - F 3.0 * s2);
          loopc4 = cT * ( F 9.0 * s2 - F 3.0 * c2);
          loopc5 = loopc1 * sT * ( F 2.0 * s2 - F 7.0 * c2);
          loopc6 = loopc1 * cT * (         c2 - F 8.0 * s2);
          for ( ; i < nt + 1; i++) {
            fi = F (i * i);  
            *ppw    = loopc2 * fi / sqrt( F 1.0 - loopc1 / fi); 
            *ppu++  = *ppw * (loopc3 + loopc5 / fi); 
            *ppw++ *=         loopc4 + loopc6 / fi;
          }

          for (i = 2; i <= btime; i++) {
            *psu++ = F 0.0;
            *psw++ = F 0.0;
          }

          a1 = F i * dt / r;
          a1 *= a1 * a1 * cons;
          f1 = ( F i + F 0.5) * srcsvel * dt / r;
          f2 = sqrt(f1 * f1 - F 1.0);
          *psw    = a1 * f2 / dt;
          *psu++  = F -3.0 * c2 * sT * *psw;
          *psw++ *= F -3.0 * s2 * cT;
          i++;

          loopc1 = r / (dt * srcsvel);
          loopc1 *= loopc1;
          loopc2 = cons * dt * dt / (r * r * r);
          loopc3 = sT * ( F 3.0 * s2 - F 9.0 * c2);
          loopc4 = cT * ( F 3.0 * c2 - F 9.0 * s2);
          loopc5 = loopc1 * sT * ( F 6.0 * c2 - F 3.0 * s2);
          loopc6 = loopc1 * cT * ( F 6.0 * s2 - F 3.0 * c2);
          for ( ; i < nt + 1; i++) {
            fi = F (i * i); 
            *psw    = loopc2 * fi / sqrt( F 1.0 - loopc1 / fi); 
            *psu++  = *psw * (loopc3 + loopc5 / fi); 
            *psw++ *=         loopc4 + loopc6 / fi;
          }
        break;
        }

        /* REDUCE TO U AND W COMPONENTS */
        ppu = tspu;
        ppw = tspw;
        psu = tssu;
        psw = tssw;
        for (i = 0; i < nt - 1; i++) {
           *ppu++ += *psu++;
           *ppw++ += *psw++;
        }

        /* CONVOLVE BY TRIANGLE OR GAUSSIAN */
        con(tsu,xx,nt-1,limit,y);
        con(tsw,xx,nt-1,limit,y);
        write(sf,tsu,4*(nt-1));
        write(sf,tsw,4*(nt-1));
      }
      ps = px; px =  pz; pz = ps;
      ss = sx; sx = -sz; sz = ss;
    }
    source++;
  }
  close(sf);
  /* MULTIPLEX WITH SYSTEM COMMAND DEMULT */
  sprintf(line, "demult nt=%d noff=%d in=%s > %s",
    2*space, nt-1, tempfile, sourcefile);
  fprintf(stderr,"nt=%d noff=%d\n",2*space,nt-1);
  system(line);
  sprintf(line, "rm %s", tempfile);
  system(line);
}

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++;
}