/*
Metaloop.c / OpenFDTD / ch.20
*/

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

int main(void)
{
	const int mode = 0;   // 0:Balmex, 1:UnBalmex
	const double fstart = 2.0e9;
	const double fstop  = 4.0e9;
	const int fdiv = 100;
	const double w = 4.0e-3;
	const double B = 3.2e-3;
	const double er = 2.6;
	const double p = 10e-3;
	const double p0 = 4.0e-3;   // 修正
	const double dg = 1.0e-3;   // 修正
	const double Cz = 0.55e-12;
	const double Ly = 2.5e-9;
	const int M = 6;   // even
	const double GPx = 110e-3;
	const double GPy = 110e-3;
	const double d = 1e-3;   // 基本セルサイズ
	const double Z0 = 80;
	double x, y;
	char title[BUFSIZ], fout[BUFSIZ];
	char smode[][BUFSIZ] = {"Balmex", "UnBalmex"};

	sprintf(title, "Metaloop Antenna M=%d d=%gmm (%s)", M, d * 1e3, smode[mode]);
	sprintf(fout, "Metaloop_M%d_d%gmm_%s.ofd", M, d * 1e3, smode[mode]);

	// initialize

	ofd_init();

	// title

	ofd_title(title);

	// mesh

	const double x0 = -GPx / 2;
	const double x1 = +GPx / 2;
	ofd_xsection(2, x0, x1);
	ofd_xdivision(1, NINT(x1 - x0, d));

	const double y0 = -GPy / 2;
	const double y1 = +GPy / 2;
	ofd_ysection(2, y0, y1);
	ofd_ydivision(1, NINT(y1 - y0, d));

	const double z0 = -2 * d;
	const double z1 = 0;
	const double z2 = z1 + B;
	const double z3 = 20e-3;
	ofd_zsection(4, z0, z1, z2, z3);
	ofd_zdivision(3, NINT(z1 - z0, d), NINT(z2 - z1, d), NINT(z3 - z2, d));

	// material

	ofd_material(er, 0.0, 1.0, 0.0, "基板");

	// geometry

	// substrate
	ofd_geometry(2, 1, x0, x1, y0, y1, z1, z2);

	// ground
	ofd_geometry(1, 1, x0, x1, y0, y1, z0, z1);

	// metaloop
	for (int side = -1; side <= +1; side += 2) {
		// X direction
		const double yw0 = side * M / 2.0 * p;
		const double yw1 = yw0 - w / 2;
		const double yw2 = yw0 + w / 2;

		// metaloop
		x = -(M / 2.0) * p;
		for (int i = 0; i < M; i++) {
			// patch
			ofd_geometry(1, 1, x, x + p0 / 2, yw1, yw2, z2, z2);
			// C
			x += p0 / 2;
			ofd_geometry(1, 1, x, x + dg, yw0, yw0, z2, z2);
			ofd_load('X', x + 0.1 * d, yw0, z2, 'C', 2 * Cz);
			// patch
			x += dg;
			ofd_geometry(1, 1, x, x + p0, yw1, yw2, z2, z2);
			// C
			x += p0;
			ofd_geometry(1, 1, x, x + dg, yw0, yw0, z2, z2);
			ofd_load('X', x + 0.1 * d, yw0, z2, 'C', 2 * Cz);
			// patch
			x += dg;
			ofd_geometry(1, 1, x, x + p0 / 2, yw1, yw2, z2, z2);
			x += p0 / 2;
		}

		// feed or L
		x = -M / 2.0 * p;
		for (int i = 0; i <= M; i++) {
			if ((side == -1) && (i == M / 2)) {
				// feed
				if (mode == 0) {
					// Balmex
					// gap
					ofd_geometry(2, 1, x - 0.9 * d, x + 0.9 * d, yw1 - 0.1 * d, yw2 + 0.1 * d, z2, z2);
					// feed
					ofd_geometry(1, 1, x - d, x + d, yw0, yw0, z2, z2);
					ofd_feed('X', x - 0.1 * d, yw0, z2, 1, 0, Z0 / 2);
					ofd_feed('X', x + 0.1 * d, yw0, z2, 1, 0, Z0 / 2);
				}
				else {
					// Unbalmex
					ofd_geometry(1, 1, x, x, yw0, yw0, z1, z2);
					ofd_feed('Z', x, yw0, z1 + 0.1 * d, 1, 0, Z0);
				}
			}
			else {
				// L
				// via
				ofd_geometry(1, 1, x, x, yw0, yw0, z1, z2);
				// slot
				ofd_geometry(2, 1, x - 0.9 * d, x + 0.9 * d, yw0 - 0.9 * d, yw0 + 0.9 * d, z1, z1);
				// L
				ofd_geometry(1, 1, x, x, yw0 + side * d, yw0, z1, z1);
				ofd_load('Y', x, yw0 + side * 0.1 * d, z1, 'L', Ly);
			}
			x += p;
		}

		// Y direction
		const double xw0 = side * M / 2.0 * p;
		const double xw1 = xw0 - w / 2;
		const double xw2 = xw0 + w / 2;

		// metaloop
		y = -(M / 2.0) * p;
		for (int j = 0; j < M; j++) {
			// patch
			ofd_geometry(1, 1, xw1, xw2, y, y + p0 / 2, z2, z2);
			// C
			y += p0 / 2;
			ofd_geometry(1, 1, xw0, xw0, y, y + dg, z2, z2);
			ofd_load('Y', xw0, y + 0.1 * d, z2, 'C', 2 * Cz);
			// patch
			y += dg;
			ofd_geometry(1, 1, xw1, xw2, y, y + p0, z2, z2);
			// C
			y += p0;
			ofd_geometry(1, 1, xw0, xw0, y, y + dg, z2, z2);
			ofd_load('Y', xw0, y + 0.1 * d, z2, 'C', 2 * Cz);
			// patch
			y += dg;
			ofd_geometry(1, 1, xw1, xw2, y, y + p0 / 2, z2, z2);
			y += p0 / 2;
		}

		// L
		y = -(M / 2.0 - 1) * p;
		for (int i = 0; i < M - 1; i++) {
			// via
			ofd_geometry(1, 1, xw0, xw0, y, y, z1, z2);
			// slot
			ofd_geometry(2, 1, xw0 - 0.9 * d, xw0 + 0.9 * d, y - 0.9 * d, y + 0.9 * d, z1, z1);
			// L
			ofd_geometry(1, 1, xw0 + side * d, xw0, y, y, z1, z1);
			ofd_load('X', xw0 + side * 0.1 * d, y, z1, 'L', Ly);
			y += p;
		}
	}

	// internal resistance
	ofd_rfeed(10);

	// ABC

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

	// frequency

	ofd_frequency1(fstart, fstop, fdiv);

	ofd_frequency2(2e9, 4e9, 40);

	// solver

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

	// iteration

	ofd_plotiter(1);

	// frequency

	//ofd_plotsmith();
	ofd_plotzin(1, 0, 0, 0);
	//ofd_plotyin(1, 0, 0, 0);
	ofd_plotref(1, 0, 0, 0);

	// far-1d

	ofd_plotfar1d('Y', 180, 0);
	ofd_far1dcomponent(1, 0, 0);
	ofd_far1dscale(-20, +15, 7);

	// window

	ofd_window3d(400, 400, 12, 60, 30);

	// output

	ofd_outdata(fout);

	return 0;
}