/*
Fan.c / OpenFDTD / ch.6
*/

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

static double radius(int shape, double z, double rp, double zp, double zc, double rc)
{
	double r;

	if (z < zp) {
		if (shape == 0) {
			r = rp * (1 - fabs((z - zp) / zp));  // straight
		}
		else {
			r = rp * (1 - (z - zp) * (z - zp) / (zp * zp));  // parabola
		}
	}
	else {
		//r = sqrt((rc * rc) - (z - zc) * (z - zc));  // circle
		r = rp * sqrt(1 - ((z - zp) * (z - zp)) / ((zc + rc - zp) * (zc + rc - zp)));  // ellipse
	}

	return r;
}

int main(void)
{
	const int shape = 1;  // 0=triangle, 1=Fan, 2=X-Fan, 3=X-Fan-W
	const double fstart =  1.0e9;
	const double fstop  = 11.0e9;
	const int fdiv = 100;
	const double L4 = 75e-3;
	const double d = 1e-3;           // cell size
	const double delta = L4 / 400;
	const double xp = d * NINT(100 * delta, d);
	const double zp = d * NINT(176 * delta, d);
	const double rc = d * NINT(116 * delta, d);
	const double zf = 1 * d;
	const double LGP = d * NINT(3 * xp, d);
	//printf("%f %f %f\n", xp, zp, rc);

	const char title[][BUFSIZ] = {
		"triangle monopole antenna with a round edge",
		"Fan-shaped antenna",
		"Cross Fan-shaped antenna",
		"Cross Fan-shaped antenna surrounded by a wire"
	};
	const char fout[][BUFSIZ] = {
		"Fan_triangle.ofd",
		"Fan.ofd",
		"Fan_X.ofd",
		"Fan_X_W.ofd"
	};

	// initialize

	ofd_init();

	// title

	ofd_title(title[shape]);

	// mesh

	const double x1 = -d * NINT(-LGP / 2, d);
	const double x2 = -x1;
	int xdiv = NINT(x2 - x1, d);
	if (xdiv % 2 == 1) xdiv++;  // even
	ofd_xsection(2, x1, x2);
	ofd_xdivision(1, xdiv);

	const double y1 = -d * NINT(-LGP / 2, d);
	const double y2 = -y1;
	int ydiv = NINT(y2 - y1, d);
	if (ydiv % 2 == 1) ydiv++;  // even
	ofd_ysection(2, y1, y2);
	ofd_ydivision(1, ydiv);

	const double z1 = -2 * d;
	const double z2 = 0;
	const double z3 = d * NINT(1.6 * zp, d);
	ofd_zsection(3, z1, z2, z3);
	ofd_zdivision(2, NINT(z2 - z1, d), NINT(z3 - z2, d));

	// geometry

	// ground
	ofd_geometry(1, 1, x1, x2, y1, y2, z1, z2);

	// feed
	ofd_geometry(1, 1, 0, 0, 0, 0, z2, z2 + zf);

	// fan
	const double zc = d * NINT(zp - sqrt((rc * rc) - (xp * xp)), d);
	const int nz = NINT(zc + rc - zf, d);
	for (int iz = 0; iz < nz; iz++) {
		const double za = zf + (iz + 0) * d;
		const double zb = zf + (iz + 1) * d;
		const double r = radius(shape, (za + zb) / 2, xp, zp, zc, rc);
		//r = d * NINT(r, d);
		ofd_geometry(1, 1, -r, +r, 0, 0, za, zb);
		if (shape > 1) {
			ofd_geometry(1, 1, 0, 0, -r, +r, za, zb);
		}
	}

	// wire
	if (shape == 3) {
		const double a = d * (int)(xp / 2 / d);
		const double b = 2 * a;
		ofd_geometry(1, 1, +b, +b, -a, +a, zp, zp);
		ofd_geometry(1, 1, -b, -b, -a, +a, zp, zp);
		ofd_geometry(1, 1, -a, +a, +b, +b, zp, zp);
		ofd_geometry(1, 1, -a, +a, -b, -b, zp, zp);
		ofd_geometry(1, 1, +b, +a, +a, +a, zp, zp);
		ofd_geometry(1, 1, +a, +a, +a, +b, zp, zp);
		ofd_geometry(1, 1, -a, -a, +b, +a, zp, zp);
		ofd_geometry(1, 1, -a, -b, +a, +a, zp, zp);
		ofd_geometry(1, 1, -b, -a, -a, -a, zp, zp);
		ofd_geometry(1, 1, -a, -a, -a, -b, zp, zp);
		ofd_geometry(1, 1, +a, +a, -b, -a, zp, zp);
		ofd_geometry(1, 1, +a, +b, -a, -a, zp, zp);
	}

	// feed

	ofd_feed('Z', 0, 0, z2 + 0.1 * d, 1, 0, 50);
	//ofd_rfeed(10);

	// frequency

	ofd_frequency1(fstart, fstop, fdiv);

	const double f0  = 3e9;
	ofd_frequency2(f0, f0, 0);

	// ABC

	ofd_pml(6, 3, 1e-6);

	// solver

	ofd_solver(20000, 200, 1e-3);

	// plot iteration

	ofd_plotiter(1);

	// frequency char.

	//ofd_plotsmith();
	ofd_plotzin(1, 0, 0, 0);
	//ofd_plotyin(1, 0, 0, 0);
	ofd_plotref(2, -40, 0, 4);
	ofd_freqdiv(10);

	// output

	ofd_outdata(fout[shape]);

	return 0;
}