/*
Fan.c / OpenMOM/ ch.6
*/

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

int main(void)
{
	const int shape = 2;  // 0=triangle, 1=Fan, 2=X-Fan, 3=X-Fan-W
	const double fstart =  1.0e9;
	const double fstop  = 11.0e9;
	const int fdiv = 50;
	const double d = 75e-3 / 400;
	const double xp = 100 * d;
	const double zp = 176 * d;
	const double rc = 116 * d;
	double dz = 3e-3;

	const double pi = 4 * atan(1);
	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.omm",
		"Fan.omm",
		"Fan_X.omm",
		"Fan_X_W.omm"
	};

	// initialize

	omm_init();

	// title

	omm_title(title[shape]);

	// geometry

	const int nz = NINT(zp, dz);
	dz = zp / nz;
	const double dx = xp / (nz - 1);

	// feed
	omm_geometry_zline(0, dz, 0, 0, 1);
	omm_feed(1, 0);
	omm_radius(dz / 2);  // 入力インピーダンスに大きく影響する

	// center line
	omm_geometry_wire(1, 0, 0, 0, 0, dz, zp, nz - 1);

	// x-lines
	const double a = (shape == 0) ? (zp - dz) / xp : xp / (dz - zp) / (dz - zp);
	for (int i = 1; i < nz; i++) {
		const double z = (i + 1) * dz;
		const double x = (shape == 0) ? (z - dz) / a : xp - a * (z - zp) * (z - zp);
		omm_geometry_wire(1, -x, +x, 0, 0, z, z, 2 * i);
		if (shape > 1) {
			const double y = x;
			omm_geometry_wire(1, 0, 0, -y, +y, z, z, 2 * i);
		}
	}

	// z-lines
	if (shape == 0) {
		for (int i = 1; i < nz - 1; i++) {
			const double x = i * dx;
			const double z = a * x + dz;
			omm_geometry_wire(1, +x, +x, 0, 0, z, zp, nz - i - 1);
			omm_geometry_wire(1, -x, -x, 0, 0, z, zp, nz - i - 1);
		}
		// edge
		omm_geometry_wire(1, 0, +xp, 0, 0, dz, zp, nz - 1);
		omm_geometry_wire(1, 0, -xp, 0, 0, dz, zp, nz - 1);
	}
	else {
		for (int i = 1; i < nz - 1; i++) {
			for (int j = 0; j < nz - i - 1; j++) {
				const double z0 = zp - (j + 0) * dz;
				const double z1 = zp - (j + 1) * dz;
				const double x0 = (xp - a * (z0 - zp) * (z0 - zp)) * i / (nz - j - 1);
				const double x1 = (xp - a * (z1 - zp) * (z1 - zp)) * i / (nz - j - 2);
				omm_geometry_wire(1, +x0, +x1, 0, 0, z0, z1, 1);
				omm_geometry_wire(1, -x0, -x1, 0, 0, z0, z1, 1);
				if (shape > 1) {
					const double y0 = x0;
					const double y1 = x1;
					omm_geometry_wire(1, 0, 0, +y0, +y1, z0, z1, 1);
					omm_geometry_wire(1, 0, 0, -y0, -y1, z0, z1, 1);
				}
			}
		}
		// edge
		for (int j = 0; j < nz - 1; j++) {
			const double z0 = zp - (j + 0) * dz;
			const double z1 = zp - (j + 1) * dz;
			const double x0 = xp - a * (z0 - zp) * (z0 - zp);
			const double x1 = xp - a * (z1 - zp) * (z1 - zp);
			omm_geometry_wire(1, +x0, +x1, 0, 0, z0, z1, 1);
			omm_geometry_wire(1, -x0, -x1, 0, 0, z0, z1, 1);
			if (shape > 1) {
				const double y0 = x0;
				const double y1 = x1;
				omm_geometry_wire(1, 0, 0, +y0, +y1, z0, z1, 1);
				omm_geometry_wire(1, 0, 0, -y0, -y1, z0, z1, 1);
			}
		}
	}

	// cap
	const double zc = zp - sqrt((rc * rc) - (xp * xp));
	omm_geometry_wire(1, 0, 0, 0, 0, zp, zc + rc, nz - 1);
	for (int i = 0; i < nz - 1; i++) {
		const double d0 = (double)(i + 0) / (nz - 1);
		const double d1 = (double)(i + 1) / (nz - 1);
		const double xe0 = xp             * d0;
		const double ze0 = (zc + rc - zp) * d0;
		const double xe1 = xp             * d1;
		const double ze1 = (zc + rc - zp) * d1;
		const double da0 = (pi / 2) / (i + 0);
		const double da1 = (pi / 2) / (i + 1);
		for (int j = 0; j <= i; j++) {
			const double arg0 = (j + 0) * da1;
			const double arg1 = (j + 1) * da1;
			omm_geometry_wire(1, -xe1 * cos(arg0), -xe1 * cos(arg1), 0, 0, zp + ze1 * sin(arg0), zp + ze1 * sin(arg1), 1);
			omm_geometry_wire(1, +xe1 * cos(arg0), +xe1 * cos(arg1), 0, 0, zp + ze1 * sin(arg0), zp + ze1 * sin(arg1), 1);
			if (shape > 1) {
				const double ye1 = xe1;
				omm_geometry_wire(1, 0, 0, -ye1 * cos(arg0), -ye1 * cos(arg1), zp + ze1 * sin(arg0), zp + ze1 * sin(arg1), 1);
				omm_geometry_wire(1, 0, 0, +ye1 * cos(arg0), +ye1 * cos(arg1), zp + ze1 * sin(arg0), zp + ze1 * sin(arg1), 1);
			}
		}
		for (int j = 1; j <= i; j++) {
			const double arg0 = (j - 1) * da0;
			const double arg1 = (j + 0) * da1;
			omm_geometry_wire(1, -xe0 * cos(arg0), -xe1 * cos(arg1), 0, 0, zp + ze0 * sin(arg0), zp + ze1 * sin(arg1), 1);
			omm_geometry_wire(1, +xe0 * cos(arg0), +xe1 * cos(arg1), 0, 0, zp + ze0 * sin(arg0), zp + ze1 * sin(arg1), 1);
			if (shape > 1) {
				const double ye0 = xe0;
				const double ye1 = xe1;
				omm_geometry_wire(1, 0, 0, -ye0 * cos(arg0), -ye1 * cos(arg1), zp + ze0 * sin(arg0), zp + ze1 * sin(arg1), 1);
				omm_geometry_wire(1, 0, 0, +ye0 * cos(arg0), +ye1 * cos(arg1), zp + ze0 * sin(arg0), zp + ze1 * sin(arg1), 1);
			}
		}
	}

	// stop wire
	if (shape == 3) {
		omm_geometry_wire(2, xp, xp, 0, 360, zp, zp, 4 * NINT(xp * pi / 2, dz));
	}

	// ground

	omm_ground();
/*
	// finite ground plane
	const int div = 2 * 10; 
	const double hw = (div / 2) * dz;
	omm_geometry_zrect(0, -hw, +hw, -hw, +hw, div, div);
*/
	// radius
	omm_radiusall(1, 0.1 * dz);

	// frequency

	omm_frequency(fstart, fstop, fdiv);

	// frequency char.

	omm_plotfrequency(1, 1, 1, 1);

	// far1d field

	//omm_plotfar1d('X', 180, 0);

	//omm_far1dstyle(0);
	//omm_far1dcomponent(1, 0, 0);
	//omm_far1ddb(1);
	//omm_far1dscale(-30, +10, 4);

	// output

	omm_outdata(fout[shape]);

	return 0;
}