/* esubs.c --- puts analytical source into emodel */
#include <stdio.h>
#include <math.h>
#define PHI theta
#define PI 3.1415926

struct disp
{
	float u; 
	float w;
};

mkring(p1,sourcew,s0,zs,nx,ps,px)
int sourcew, s0, zs;
int nx;
float *ps;
struct disp *px, *p1;
{
	int i;
	float *pps, theta;

	pps = ps;
	p1 = p1 + s0 - sourcew + (zs - sourcew) * nx;
	for (i=0; i<= sourcew; i++)
	{
		theta = -atan((sourcew-i) / (float)sourcew);
		px->u = p1->u - *pps * sin(PHI);
		px->w = p1->w + *pps * cos(PHI);
		p1++;
		pps++;
		px++;
	}
	pps = pps - 2;
	for (i=sourcew+1; i<2*sourcew; i++)
	{
		theta = atan((i-sourcew) / (float)sourcew); 
		px->u = p1->u - *pps * sin(PHI);
		px->w = p1->w + *pps * cos(PHI);
		p1++;
		pps--;
		px++;
	}
	pps = ps;
	for (i=2*sourcew; i<=3*sourcew; i++)
	{
		theta = PI/2. - atan((3*sourcew-i) / (float)sourcew);
		px->u = p1->u - *pps * sin(PHI);
		px->w = p1->w + *pps * cos(PHI);
		p1 += nx;
		pps++;
		px++;
	}
	pps = pps - 2;
	for (i=3*sourcew+1; i<4*sourcew; i++)
	{
		theta = PI/2. + atan((i-3*sourcew) / (float)sourcew); 
		px->u = p1->u - *pps * sin(PHI);
		px->w = p1->w + *pps * cos(PHI);
		p1 += nx;
		pps--;
		px++;
	}
	pps = ps;
	for (i=4*sourcew; i<=5*sourcew; i++)
	{
		theta = PI - atan((5*sourcew-i) / (float)sourcew);
		px->u = p1->u - *pps * sin(PHI);
		px->w = p1->w + *pps * cos(PHI);
		p1--;
		pps++;
		px++;
	}
	pps = pps - 2;
	for (i=5*sourcew+1; i<6*sourcew; i++)
	{
		theta = PI + atan((i-5*sourcew) / (float)sourcew); 
		px->u = p1->u - *pps * sin(PHI);
		px->w = p1->w + *pps * cos(PHI);
		p1--;
		pps--;
		px++;
	}
	pps = ps;
	for (i=6*sourcew; i<=7*sourcew; i++)
	{
		theta = 3.*PI/2. - atan((7*sourcew-i) / (float)sourcew);
		px->u = p1->u - *pps * sin(PHI);
		px->w = p1->w + *pps * cos(PHI);
		p1 -= nx;
		pps++;
		px++;
	}
	pps = pps - 2;
	for (i=7*sourcew+1; i<8*sourcew; i++)
	{
		theta = 3.*PI/2. + atan((i-7*sourcew) / (float)sourcew); 
		px->u = p1->u - *pps * sin(PHI);
		px->w = p1->w + *pps * cos(PHI);
		p1 -= nx;
		pps--;
		px++;
	}
}

mkring2(p1,sourcew,s0,zs,nx,px)
int sourcew, s0, zs;
int nx;
struct disp *p1, *px;
{
	int i;

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

bkring(p1,sourcew,s0,zs,nx,pfg)
int sourcew, s0, zs;
int nx;
struct disp *pfg, *p1 ;
{
	int i;

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

addsource(p1,sourcew,s0,zs,nx,ps)
int sourcew, s0, zs;
int nx;
float *ps;
struct disp *p1;
{
	int i;
	float *pps, theta;

	pps = ps;
	p1 = p1 + s0 - sourcew +(zs - sourcew) * nx;
	for (i=0; i<=sourcew; i++)
	{
		theta = -atan((sourcew-i) / (float)sourcew);
		p1->u += *pps * sin(PHI);
		p1->w -= *pps * cos(PHI);
		p1++;
		pps++;
	}
	pps = pps - 2;
	for (i=sourcew+1; i<2*sourcew; i++)
	{
		theta = atan((i-sourcew) / (float)sourcew); 
		p1->u += *pps * sin(PHI);
		p1->w -= *pps * cos(PHI);
		p1++;
		pps--;
	}
	pps = ps;
	for (i=2*sourcew; i<=3*sourcew; i++)
	{
		theta = PI/2. - atan((3*sourcew-i) / (float)sourcew);
		p1->u += *pps * sin(PHI);
		p1->w -= *pps * cos(PHI);
		p1 += nx;
		pps++;
	}
	pps = pps - 2;
	for (i=3*sourcew+1; i<4*sourcew; i++)
	{
		theta = PI/2. + atan((i-3*sourcew) / (float)sourcew); 
		p1->u += *pps * sin(PHI);
		p1->w -= *pps * cos(PHI);
		p1 += nx;
		pps--;
	}
	pps = ps;
	for (i=4*sourcew; i<=5*sourcew; i++)
	{
		theta = PI - atan((5*sourcew-i) / (float)sourcew);
		p1->u += *pps * sin(PHI);
		p1->w -= *pps * cos(PHI);
		p1--;
		pps++;
	}
	pps = pps - 2;
	for (i=5*sourcew+1; i<6*sourcew; i++)
	{
		theta = PI + atan((i-5*sourcew) / (float)sourcew); 
		p1->u += *pps * sin(PHI);
		p1->w -= *pps * cos(PHI);
		p1--;
		pps--;
	}
	pps = ps;
	for (i=6*sourcew; i<=7*sourcew; i++)
	{
		theta = 3.*PI/2. - atan((7*sourcew-i) / (float)sourcew);
		p1->u += *pps * sin(PHI);
		p1->w -= *pps * cos(PHI);
		p1 -= nx;
		pps++;
	}
	pps = pps - 2;
	for (i=7*sourcew+1; i<8*sourcew; i++)
	{
		theta = 3.*PI/2. + atan((i-7*sourcew) / (float)sourcew); 
		p1->u += *pps * sin(PHI);
		p1->w -= *pps * cos(PHI);
		p1 -= nx;
		pps--;
	}
}

savit(p3,source,s0,zs,nx,px)
struct disp *p3, *px;
int source, s0, zs, nx;
{
	struct disp *pp3;
	int iz, ix;

	for (iz=zs-source-1; iz<= zs+source+1; iz++)
	{
		pp3 = p3 + s0 - source - 1 + nx * iz;
		for (ix=s0-source-1; ix<=s0+source+1; ix++)
		{
			px->u = pp3->u;
			px->w = pp3->w;
			pp3++;
			px++;
		}
	}
}

recall(p3,source,s0,zs,nx,px)
struct disp *p3, *px;
int source, s0, zs, nx;
{
	struct disp *pp3;
	int iz, ix;

	for (iz=zs-source-1; iz<= zs+source+1; iz++)
	{
		pp3 = p3 + s0 - source - 1 + nx * iz;
		for(ix=s0-source-1; ix<=s0+source+1; ix++) 
		{
			pp3->u = px->u;
			pp3->w = px->w;
			pp3++;
			px++;
		}
	}
}