/*
CP_GAA.c / OpenMOM / ch.15
*/

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

static void loop(int bowtie, double x0, double y0, double z0, double lx, double r, double ds, double d, double dy)
{
	// loop
	const double pi = 4 * atan(1);
	const double dtor = pi / 180;
	const double a = atan(d / 2 / r);
	const double l = 2 * (pi - a) * r;
	int div = NINT(l, d);
	if (div % 2 == 1) div++;  // even
	omm_geometry_wire(2, r, r, 90 + a / dtor, 450 - a / dtor, 0, 0, div);
	const double yc = y0 - ds - r;
	omm_offset(x0, yc, z0);

	// 摂動
	if (bowtie) {
		const double p = 3e-3;
		omm_geometry_yline(yc - r, yc - r + p, z0, x0, NINT(p, d));
	}

	// Y-line
	const double dx = r * sin(a);
	const double y1 = yc + r * cos(a);
	omm_geometry_yline(y0, y1, z0, x0 + dx, NINT(y1 - y0, d));
	omm_geometry_yline(y0, y1, z0, x0 - dx, NINT(y1 - y0, d));

	// X-line
	int div2 = NINT(lx / 2 - dx, d);
	if (bowtie) {
		omm_geometry_plane(1, x0 + dx, x0 + lx / 2, x0 + lx / 2, x0 + dx,
		                      y0, y0 + dy, y0 - dy, y0, z0, z0, z0, z0, div2, 1);
		omm_geometry_plane(1, x0 - dx, x0 - lx / 2, x0 - lx / 2, x0 - dx,
		                      y0, y0 + dy, y0 - dy, y0, z0, z0, z0, z0, div2, 1);
	}
	else {
		omm_geometry_xline(x0 + dx, x0 + lx / 2, y0, z0, div2);
		omm_geometry_xline(x0 - dx, x0 - lx / 2, y0, z0, div2);
	}
}

int main(void)
{
	const double fstart = 5.0e9;
	const double fstop  = 7.5e9;
	const int fdiv = 50;
	const int N = 11;  // odd
	const double Lx = 25e-3;
	const double Ly = 50e-3;
	const double H = 2.5e-3;
	const double w1 = 2.8e-3;
	const double w2 = 0.8e-3;
	const double ds = 2.5e-3;
	const double R0 = 50;
	const double d = 2.5e-3;   // 基本要素長
	const double r1 = w1 / 4;
	const double r2 = w2 / 4;

	// variables
	const double rlp = 9.5e-3;

	// initialize

	omm_init();

	// title

	omm_title("Circularly Polarized Grid Array Antenna");

	// geometry

	// ground
	omm_ground();

	// feed
	const double x0 = - (N / 2.0) * Lx;
	omm_geometry_zline(0, H, + x0, 0, NINT(H, d));
	omm_feed(1, 0);
	omm_radius(r1);

	// feed
	omm_geometry_yline(- Ly / 2 + w1 / 2, + Ly / 2 - w1 / 2, H, + x0, NINT(Ly - w1, d));
	omm_radius(r1);

	// load
	omm_geometry_yline(- Ly / 2 + w1 / 2, + Ly / 2 - w1 / 2, H, - x0, NINT(Ly - w1, d));
	omm_radius(r1);

	// load
	omm_geometry_zline(0, H, - x0, 0, NINT(H, d));
	omm_load(1, R0);
	omm_radius(r1);

	// Y-line
	for (int i = 1; i < N; i++) {
		const double x = x0 + (i * Lx);
		omm_geometry_yline(- Ly,              - Ly / 2 - w1 / 2, H, x, NINT(Ly / 2 - w1 / 2, d));
		omm_geometry_yline(- Ly / 2 + w1 / 2,          - w1 / 2, H, x, NINT(Ly / 2 - w1,     d));
		omm_geometry_yline(+ Ly / 2 - w1 / 2,          + w1 / 2, H, x, NINT(Ly / 2 - w1,     d));
		omm_geometry_yline(+ Ly,              + Ly / 2 + w1 / 2, H, x, NINT(Ly / 2 - w1 / 2, d));
	}

	// loop
	for (int i = 0; i < N; i += 2) {
		const double x = x0 + (i + 0.5) * Lx;
		for (int side = -1; side <= 1; side += 2) {
			const double y = side * (Ly / 2);
			loop(1, x, y, H, Lx, rlp, ds, d, w1 / 2);
		}
	}
	for (int i = 1; i < N; i += 2) {
		const double x = x0 + (i + 0.5) * Lx;
		for (int side = -1; side <= 1; side++) {
			const double y = side * Ly;
			loop((side == 0), x, y, H, Lx, rlp, ds, d, w1 / 2);
		}
	}

	// radius
	omm_radiusall(1, r2);

	// frequency

	omm_frequency(fstart, fstop, fdiv);

	// frequency char.

	//omm_plotsmith();
	omm_plotzin(1, 0, 0, 0);
	omm_plotyin(1, 0, 0, 0);
	omm_plotref(1, 0, 0, 0);
	//omm_freqdiv(5);

	// far1d field

	//omm_plotfar1d('Y', 360, 0);
	//omm_far1dcomponent(0, 1, 0);
	//omm_far1dscale(-10, +25, 7);

	// window

	omm_window3d(400, 400, 12, 0, -90);

	// output

	omm_outdata("CP_GAA.omm");

	return 0;
}