#include  <stdio.h>
#include  <math.h>

#define PHI theta

extern int style;

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

	pps =ps;
	pi=3.1415926;
	p1= p1 + s0- sourcew +(zs-sourcew)*nx;
	for (i=0; i<= sourcew; i++)
		{
		theta = - atan((sourcew-i)/(float)sourcew);
		if(style==0) *px = *p1 - *pps;
		if(style==1) *px = *p1 - *pps * sin(PHI);
		if(style==2) *px = *p1 - *pps * cos(PHI);
		p1 ++ ;
		pps ++ ;
		px ++ ;
		}
	pps= pps-2;
	for (i=sourcew+1; i<2*sourcew; i++)
		{
		theta = atan((i-sourcew)/(float)sourcew); 
		if(style==0) *px = *p1 - *pps;
		if(style==1) *px = *p1 - *pps * sin(PHI);
		if(style==2) *px = *p1 - *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);
		if(style==0) *px = *p1 - *pps;
		if(style==1) *px = *p1 - *pps * sin(PHI);
		if(style==2) *px = *p1 - *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); 
		if(style==0) *px = *p1 - *pps;
		if(style==1) *px = *p1 - *pps * sin(PHI);
		if(style==2) *px = *p1 - *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);
		if(style==0) *px = *p1 - *pps;
		if(style==1) *px = *p1 - *pps * sin(PHI);
		if(style==2) *px = *p1 - *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); 
		if(style==0) *px = *p1 - *pps;
		if(style==1) *px = *p1 - *pps * sin(PHI);
		if(style==2) *px = *p1 - *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);
		if(style==0) *px = *p1 - *pps;
		if(style==1) *px = *p1 - *pps * sin(PHI);
		if(style==2) *px = *p1 - *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); 
		if(style==0) *px = *p1 - *pps;
		if(style==1) *px = *p1 - *pps * sin(PHI);
		if(style==2) *px = *p1 - *pps * cos(PHI);
		p1 -= nx;
		pps --;
		px ++;
		}
	}

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

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

bkring(p1,sourcew,s0,zs,nx,pfg)
int sourcew,s0,zs;
float *pfg;
int nx;
float *p1 ;
	{
	int i;
	p1= p1+s0-sourcew + (zs-sourcew)*nx;
	for (i=0; i<2*sourcew; i++)
		{
		*p1= *pfg;
		p1 ++;
		pfg ++;
		}
	for (i=2*sourcew; i<4*sourcew; i++)
		{
		*p1 = *pfg;
		p1 += nx;
		pfg ++;
		}
	for (i=4*sourcew; i<6*sourcew; i++)
		{
		*p1= *pfg;
		p1 --;  
		pfg ++;
		}
	for (i=6*sourcew; i<8*sourcew; i++)
		{
		*p1= *pfg;
		p1 -= nx;
		pfg ++;
		}
	}

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

	pi=3.1415926;
	pps =ps;
	p1= p1 + s0-sourcew +(zs-sourcew)*nx;
	for (i=0; i<=sourcew; i++)
		{
		theta = - atan((sourcew-i)/(float)sourcew);
		if(style==0) *p1 = *p1 + *pps;
		if(style==1) *p1 = *p1 + *pps * sin(PHI);
		if(style==2) *p1 = *p1 + *pps * cos(PHI);
		p1 ++;
		pps ++;
		}
	pps = pps-2;
	for (i=sourcew+1; i<2*sourcew; i++)
		{
		theta = atan((i-sourcew)/(float)sourcew); 
		if(style==0) *p1 = *p1 + *pps;
		if(style==1) *p1 = *p1 + *pps * sin(PHI);
		if(style==2) *p1 = *p1 + *pps * cos(PHI);
		p1 ++;
		pps --;
		}
	pps = ps;
	for (i=2*sourcew; i<=3*sourcew; i++)
		{
		theta = pi/2. - atan((3*sourcew-i)/(float)sourcew);
		if(style==0) *p1 = *p1 + *pps;
		if(style==1) *p1 = *p1 + *pps * sin(PHI);
		if(style==2) *p1 = *p1 + *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); 
		if(style==0) *p1 = *p1 + *pps;
		if(style==1) *p1 = *p1 + *pps * sin(PHI);
		if(style==2) *p1 = *p1 + *pps * cos(PHI);
		p1 += nx;
		pps --;
		}
	pps =ps;
	for (i=4*sourcew; i<=5*sourcew; i++)
		{
		theta = pi - atan((5*sourcew-i)/(float)sourcew);
		if(style==0) *p1 = *p1 + *pps;
		if(style==1) *p1 = *p1 + *pps * sin(PHI);
		if(style==2) *p1 = *p1 + *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); 
		if(style==0) *p1 = *p1 + *pps;
		if(style==1) *p1 = *p1 + *pps * sin(PHI);
		if(style==2) *p1 = *p1 + *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);
		if(style==0) *p1 = *p1 + *pps;
		if(style==1) *p1 = *p1 + *pps * sin(PHI);
		if(style==2) *p1 = *p1 + *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); 
		if(style==0) *p1 = *p1 + *pps;
		if(style==1) *p1 = *p1 + *pps * sin(PHI);
		if(style==2) *p1 = *p1 + *pps * cos(PHI);
		p1 -=nx;
		pps --;
		}
	}
savit(p3,source,s0,zs,nx,px)
float *p3, *px;
int source,s0,zs,nx;
	{
	float *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= *pp3;
			pp3 ++ ;
			px ++ ;
			}
		}
	}
recall(p3,source,s0,zs,nx,px)
float *p3, *px;
int source,s0,zs,nx;
	{
	float *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= *px;
			pp3 ++;
			px ++;
			}
		}
	}