/*
Metaloop_open.c / OpenFDTD / ch.21
*/

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

int main(void)
{
	const int slab = 0;         // 0/1 : slab
	const double fstart = 2.0e9;
	const double fstop  = 4.0e9;
	const int fdiv = 100;
	const double w = 8.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.65e-12;
	const double Ly = 1.2e-9;
	const int M = 5;  // odd
	const double GPx = 110e-3;
	const double GPy = 110e-3;
	const double d = 1e-3;   // 基本セルサイズ
	const double Z0 = 50;
	const double dslab = 58.8e-3;
	const double hslab = 21.5e-3;

	double x, y, da1, da2, db1, db2;
	char title[BUFSIZ], fout[BUFSIZ];

	sprintf(title, "Metaloop_open Antenna M=%d d=%gmm%s", M, d * 1e3, (slab ? " slab" : ""));
	sprintf(fout, "Metaloop_open_M%d_d%gmm%s.ofd", M, d * 1e3, (slab ? "_slab" : ""));

	// 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 + (slab * (dslab + hslab));
	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++) {
			da1 = da2 = db1 = db2 = 0;
			if (i == 0) {
				if (side > 0) da1 = +w / 4;  // -X端,-Yカット分
				if (side < 0) da2 = -w / 4;  // -X端,+Yカット分
			}
			else if (i == M - 1) {
				if (side > 0) db1 = +w / 4;  // +X端,-Yカット分
				if (side < 0) db2 = -w / 4;  // +X端,+Yカット分
			}
			// patch
			ofd_geometry(1, 1, x, x + p0 / 2, yw1 + da1, yw2 + da2, 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 / 2, yw1 + da1, yw2 + da2, z2, z2);
			// patch
			x += p0 / 2;
			ofd_geometry(1, 1, x, x + p0 / 2, yw1 + db1, yw2 + db2, 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 / 2, yw1 + db1, yw2 + db2, z2, z2);
			x += p0 / 2;
		}

		// L
		x = -(M - 1) / 2.0 * p;
		for (int i = 0; i < M; i++) {
			// 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++) {
			if ((side == -1) || (j != (M - 1) / 2)) {
				da1 = da2 = db1 = db2 = 0;
				if (j == 0) {
					if (side > 0) da1 = +w / 4;  // -Y端,-Xカット分
					if (side < 0) da2 = -w / 4;  // -Y端,+Xカット分
				}
				else if (j == M - 1) {
					if (side > 0) db1 = +w / 4;  // +Y端,-Xカット分
					if (side < 0) db2 = -w / 4;  // +Y端,+Xカット分
				}
				// patch
				ofd_geometry(1, 1, xw1 + da1, xw2 + da2, 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 + da1, xw2 + da2, y, y + p0 / 2, z2, z2);
				// patch
				y += p0 / 2;
				ofd_geometry(1, 1, xw1 + db1, xw2 + db2, 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 + db1, xw2 + db2, y, y + p0 / 2, z2, z2);
				y += p0 / 2;
			}
			else {
				y += p;
			}
		}

		// L
		y = -(M - 1) / 2.0 * p;
		for (int j = 0; j < M; j++) {
			if ((side == -1) || (j != (M - 1) / 2)) {
				// 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;
		}
	}

	// feed
	x = (M / 2.0) * p;
	y = p / 2;
	ofd_geometry(1, 1, x, x, y, y, z1, z2);  // wire
	//const double dh = B / NINT(B, d);
	//ofd_geometry(1, 1, x, x + d, y, y + d, z1 + dh, z2);
	//ofd_geometry(1, 1, x, x, y, y, z1, z1 + dh);
	ofd_feed('Z', x, y, z1 + 0.1 * d, 1, 0, Z0);
	ofd_rfeed(10);

	// load
	x = (M / 2.0) * p;
	y = -p / 2;
	ofd_geometry(1, 1, x, x, y, y, z1, z2);
	ofd_load('Z', x, y, z1 + 0.1 * d, 'R', Z0);

	// slab
	if (slab) {
		const double zs1 = z2 + 58.8e-3;
		const double zs2 = zs1 + 21.5e-3;
		ofd_geometry(2, 1, x0, x1, y0, y1, zs1, zs2);
	}

	// 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(0, 0, 1);
	ofd_far1dscale(-20, +10, 6);

	// window

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

	// output

	ofd_outdata(fout);

	return 0;
}