/* eqsubs.c --- for earthquakes with tdesource */
#include <stdio.h>
#include <math.h>
#define Pu ps[0]
#define Pw ps[1]

struct disp
{
	float u; 
	float w;
};

mkeqring(p1,sw,s0,zs,nx,ps,px)
int sw, s0, zs;
int nx;
float *ps;
struct disp *px, *p1;
{
	int i;

	p1 = p1 + s0 - sw + (zs - sw) * nx;
	for (i=0; i< 2*sw; i++)
	{
		px->u = p1->u - Pu;
		px->w = p1->w - Pw;
		p1++;
		ps += 2;
		px++;
	}
	for (i=2*sw; i<4*sw; i++)
	{
		px->u = p1->u - Pu;
		px->w = p1->w - Pw;
		p1 += nx;
		ps += 2;
		px++;
	}
	for (i=4*sw; i<6*sw; i++)
	{
		px->u = p1->u - Pu;
		px->w = p1->w - Pw;
		p1--;
		ps += 2;
		px++;
	}
	for (i=6*sw; i<8*sw; i++)
	{
		px->u = p1->u - Pu;
		px->w = p1->w - Pw;
		p1 -= nx;
		ps += 2;
		px++;
	}
}

addeqsource(p1,sw,s0,zs,nx,ps)
int sw, s0, zs, nx;
float *ps;
struct disp *p1;
{
	int i;

	p1 = p1 + s0 - sw + (zs - sw) * nx;
	for (i=0; i<2*sw; i++)
	{
		p1->u += Pu;
		p1->w += Pw;
		p1++;
		ps += 2;
	}
	for (i=2*sw; i<4*sw; i++)
	{
		p1->u += Pu;
		p1->w += Pw;
		p1 += nx;
		ps += 2;
	}
	for (i=4*sw; i<6*sw; i++)
	{
		p1->u += Pu;
		p1->w += Pw;
		p1--;
		ps += 2;
	}
	for (i=6*sw; i<8*sw; i++)
	{
		p1->u += Pu;
		p1->w += Pw;
		p1 -= nx;
		ps += 2;
	}
}