/*
Metaline2.c / OpenFDTD / ch.19
*/

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

int main(void)
{
	const double fstart = 2.0e9;
	const double fstop  = 4.0e9;
	const int fdiv = 100;
	const double ds = 4.0e-3;
	const double w = 8.8e-3;
	const double B = 3.2e-3;
	const double er = 2.6;
	const double p = 10e-3;
	const double p0 = 4e-3;
	const double dg = 1e-3;
	const double rvia = 0.5e-3;
	const double hvia0 = 2.2e-3;
	const double Cz = 0.51e-12;
	const double Ly = 1.28e-9;
	const double RB = 50;
	const int M = 10;
	const double GPx = 120e-3;
	const double GPy = 50e-3;
	const double d = 0.5e-3;   // 基本セルサイズ
	double x;
	char title[BUFSIZ], fout[BUFSIZ];

	sprintf(title, "Two Metaline Antennas M=%d d=%gmm", M, d * 1e3);
	sprintf(fout, "Metaline2_M%d_d%gmm.ofd", M, d * 1e3);

	// 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 = -ds / 2 - w;
	const double y2 = -ds / 2 - w / 2;
	const double y3 = -ds / 2;
	const double y4 = -y3;
	const double y5 = -y2;
	const double y6 = -y1;
	const double y7 = -y0;
	ofd_ysection(8, y0, y1, y2, y3, y4, y5, y6, y7);
	ofd_ydivision(7, NINT(y1 - y0, d), NINT(y2 - y1, d), NINT(y3 - y2, d), NINT(y4 - y3, d), NINT(y5 - y4, d), NINT(y6 - y5, d), NINT(y7 - y6, 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, y7, z1, z2);

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

	// setup via
	//const double rvia = d * NINT(rvia0, d);  // 0 ?
	const double dh = B / NINT(B, d);  // 基板内厚さ方向セル
	const double hvia = dh * NINT(hvia0, dh);

	// metaline
	for (int side = -1; side <= +1; side += 2) {
		// Y center
		const double yc = side * (ds / 2 + w / 2);

		// feed
		x = -side * (M / 2.0 * p);
		ofd_geometry(1, 1, x, x, yc, yc, z1, z2);
		ofd_feed('Z', x, yc, z1 + 0.1 * d, side, 0, RB);
		ofd_rfeed(10);

		// metaline
		for (int i = 0; i < M; i++) {
			// patch
			x = (-M / 2.0 + i) * p;
			ofd_geometry(1, 1, x, x + p0 / 2, yc - w / 2, yc + w / 2, z2, z2);
			// C
			x += p0 / 2;
			ofd_geometry(1, 1, x, x + dg, yc, yc, z2, z2);
			ofd_load('X', x + 0.1 * d, yc, z2, 'C', 2 * Cz);
			// patch
			x += dg;
			ofd_geometry(1, 1, x, x + p0,     yc - w / 2, yc + w / 2, z2, z2);
			// via + L
			x += p0 / 2;
			if (NINT(rvia,  d / 2) <= 1) {  // via = 1 cell
				const double xv = x  - (side < 0 ? d : 0);
				const double yv = yc - (side < 0 ? d : 0);
				ofd_geometry(1, 1, xv, xv + d, yv, yv + d, z2 - hvia, z2);
			}
			else {
				ofd_geometry(1, 1, x - rvia, x + rvia, yc - rvia, yc + rvia, z2 - hvia, z2);
			}
			ofd_geometry(1, 1, x, x, yc, yc, z1, z2 - hvia);
			ofd_load('Z', x, yc, z1 + 0.1 * d, 'L', Ly);
			// C
			x += p0 / 2;
			ofd_geometry(1, 1, x, x + dg, yc, yc, z2, z2);
			ofd_load('X', x + 0.1 * d, yc, z2, 'C', 2 * Cz);
			// patch
			x += dg;
			ofd_geometry(1, 1, x, x + p0 / 2, yc - w / 2, yc + w / 2, z2, z2);
		}

		// load
		x = +side * (M / 2.0 * p);
		ofd_geometry(1, 1, x, x, yc, yc, z1, z2);
		ofd_load('Z', x, yc, 0.1 * d, 'R', RB);
	}

	// ABC

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

	// frequency

	ofd_frequency1(fstart, fstop, fdiv);

	ofd_frequency2(2.4e9, 3.6e9, 6);

	// solver

	ofd_solver(NINT(150, d), NINT(0.2, d), 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, +10, 6);

	// window

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

	// output

	ofd_outdata(fout);

	return 0;
}