/*
Metahelical.c / OpenFDTD / ch.23
*/

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

static void metaline_x(double, double, double, double, double, int, double, double, double, int, double);
static void metaline_y(double, double, double, double, double, int, double, double, double, int, double);

int main(void)
{
	const int N = 2;
	const double w = 2.0e-3;
	const double B = 2.0e-3;   // 修正
	const double er = 2.6;
	const double S = 40e-3;    // 修正
	const double p0 = 4e-3;
	const double dg = 1e-3;
	const double p = 2 * (p0 + dg);
	//const double rvia = 0.5e-3;
	//const double hvia = 0.6e-3;
	const double dh = 3e-3;     // 可変?
	const double h = (4 * N + 2) * dh;
	const double Cz = 0.53e-12;
	const double Ly = 3.83e-9;
	const double RB = 88;
	const double margin = 20e-3;
	const double GPx = S + (2 * margin);
	const double GPy = GPx;
	const double d = 1e-3;   // 基本セルサイズ
	char title[BUFSIZ], fout[BUFSIZ];

	sprintf(title, "Metahelical Antenna N=%d d=%gmm", N, d * 1e3);
	sprintf(fout, "Metahelical_N%d_d%gmm.ofd", N, 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 = +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 + h;
	const double z3 = z2 + margin;
	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, "基板");

	// Cをloadで表現すると物性値数が256を超えるので物性値番号3のgeometryで表現する
	const double eps0 = 8.854e-12;
	const double db = B / NINT(B, d);  // 基板中のセルサイズ
	const double ds = (d + db) / 2;  // 基板表面セルサイズの平均
	const double er3 = (2 * Cz) / eps0 * dg / (d * ds);
	ofd_material(er3, 0, 1, 0, "2Cz");

	// geometry

	// substrate
	ofd_geometry(2, 1, -S / 2, -S / 2 + B, -S / 2, +S / 2,     0, h);
	ofd_geometry(2, 1, +S / 2, +S / 2 - B, -S / 2, +S / 2,     0, h);
	ofd_geometry(2, 1, -S / 2, +S / 2,     -S / 2, -S / 2 + B, 0, h);
	ofd_geometry(2, 1, -S / 2, +S / 2,     +S / 2, +S / 2 - B, 0, h);

	// ground
	ofd_geometry(1, 1, -S / 2 + B, -S / 2 + B, -S / 2 + B, +S / 2 - B, 0, h);
	ofd_geometry(1, 1, +S / 2 - B, +S / 2 - B, -S / 2 + B, +S / 2 - B, 0, h);
	ofd_geometry(1, 1, -S / 2 + B, +S / 2 - B, -S / 2 + B, -S / 2 + B, 0, h);
	ofd_geometry(1, 1, -S / 2 + B, +S / 2 - B, +S / 2 - B, +S / 2 - B, 0, h);

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

	// feed
	ofd_geometry(1, 1, +S / 2, +S / 2, -S / 2, -S / 2, z1, z1 + d);
	ofd_feed('Z', +S / 2, -S / 2, z1 + 0.1 * d, 1, 0, RB);
	ofd_rfeed(10);

	// metahelical
	const int np = NINT(S, p);
	double z = z1 + d + dh + w / 2;
	for (int turn = 0; turn < N; turn++) {
		// +X
		metaline_y(+S / 2 - B, +S / 2, -S / 2, +S / 2, z, np, dg, w, dh, 3, Ly);
		z += dh;
		// +Y
		metaline_x(+S / 2, -S / 2, +S / 2 - B, +S / 2, z, np, dg, w, dh, 3, Ly);
		z += dh;
		// -X
		metaline_y(-S / 2 + B, -S / 2, +S / 2, -S / 2, z, np, dg, w, dh, 3, Ly);
		z += dh;
		// -Y
		metaline_x(-S / 2, +S / 2, -S / 2 + B, -S / 2, z, np, dg, w, dh, 3, Ly);
		z += dh;
	}

	// load
	const double xload = +S / 2 - B;
	const double yload = -S / 2;
	const double zload = z - dh;
	ofd_geometry(1, 1, xload, xload, yload, yload + B, zload, zload);
	ofd_load('Y', xload, yload + 0.1 * d, zload, 'R', RB);

	// ABC

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

	// frequency

	ofd_frequency1(2e9, 5e9, 150);

	ofd_frequency2(2e9, 5e9, 60);

	// 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(2, -25, 0, 5);
	ofd_freqdiv(6);

	// far-0d

	ofd_plotfar0d(0, 0, 2, -20, +10, 6);

	// 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;
}

/*
X方向
x1:始点のX座標
x2:終点のX座標
y1:基板底面のY座標
y2:基板表面のY座標
z0:中心のZ座標
np:周期の数
*/
static void metaline_x(double x1, double x2, double y1, double y2, double z0,
	int np, double dg, double w, double dh, int m, double Ly)
{
	dg = (x2 > x1) ? +dg : -dg;
	const double p = (x2 - x1) / np;
	const double p0 = (p - 2 * dg) / 2;
	const double xstep = (x2 - x1) ? +p0 / 2 : -p0 / 2;
	const double z1 = z0 - w / 2;
	const double z2 = z0 + w / 2;
/*
	ofd_geometry(1, 1, x1, x1 + p0 / 2, y, y, z0 - w / 2 - dh, z0 - w / 2);
	ofd_geometry(1, 1, x1, x2,          y, y, z0 - w / 2,      z0 + w / 2);
	ofd_geometry(1, 1, x2, x2 - p0 / 2, y, y, z0 + w / 2 + dh, z0 + w / 2);
*/
	// step
	ofd_geometry(1, 1, x1, x1 + xstep, y2, y2, z1, z1 - dh);

	double x = x1;
	for (int i = 0; i < np; i++) {
		// patch
		ofd_geometry(1, 1, x, x + p0 / 2, y2, y2, z1, z2);
		x += p0 / 2;

		// C
		ofd_geometry(m, 1, x, x + dg, y2, y2, z0, z0);
		x += dg;

		// patch
		ofd_geometry(1, 1, x, x + p0 / 2, y2, y2, z1, z2);
		x += p0 / 2;

		// via + L
		ofd_geometry(1, 1, x, x, y1, y2, z0, z0);
		ofd_load('Y', x, y1 + 0.1 * (y2 - y1), z0, 'L', Ly);

		// patch
		ofd_geometry(1, 1, x, x + p0 / 2, y2, y2, z1, z2);
		x += p0 / 2;

		// C
		ofd_geometry(m, 1, x, x + dg, y2, y2, z0, z0);
		x += dg;

		// patch
		ofd_geometry(1, 1, x, x + p0 / 2, y2, y2, z1, z2);
		x += p0 / 2;
	}

	// step
	ofd_geometry(1, 1, x2, x2 - xstep, y2, y2, z2, z2 + dh);
}

/*
Y方向
x:X座標
y1:始点のY座標
y2:終点のY座標
z0:中心のZ座標
np:周期の数
*/
static void metaline_y(double x1, double x2, double y1, double y2, double z0,
	int np, double dg, double w, double dh, int m, double Ly)
{
	dg = (y2 > y1) ? +dg : -dg;
	const double p = (y2 - y1) / np;
	const double p0 = (p - 2 * dg) / 2;
	const double ystep = (y2 - y1) ? +p0 / 2 : -p0 / 2;
	const double z1 = z0 - w / 2;
	const double z2 = z0 + w / 2;
/*
	ofd_geometry(1, 1, x, x, y1, y1 + p0 / 2, z0 - w / 2 - dh, z0 - w / 2);
	ofd_geometry(1, 1, x, x, y1, y2,          z0 - w / 2,      z0 + w / 2);
	ofd_geometry(1, 1, x, x, y2, y2 - p0 / 2, z0 + w / 2 + dh, z0 + w / 2);
*/
	// step
	ofd_geometry(1, 1, x2, x2, y1, y1 + ystep, z1, z1 - dh);

	double y = y1;
	for (int i = 0; i < np; i++) {
		// patch
		ofd_geometry(1, 1, x2, x2, y, y + p0 / 2, z1, z2);
		y += p0 / 2;

		// C
		ofd_geometry(m, 1, x2, x2, y, y + dg, z0, z0);
		y += dg;

		// patch
		ofd_geometry(1, 1, x2, x2, y, y + p0 / 2, z1, z2);
		y += p0 / 2;

		// via + L
		ofd_geometry(1, 1, x1, x2, y, y, z0, z0);
		ofd_load('X', x1 + 0.1 * (x2 - x1), y, z0, 'L', Ly);

		// patch
		ofd_geometry(1, 1, x2, x2, y, y + p0 / 2, z1, z2);
		y += p0 / 2;

		// C
		ofd_geometry(m, 1, x2, x2, y, y + dg, z0, z0);
		y += dg;

		// patch
		ofd_geometry(1, 1, x2, x2, y, y + p0 / 2, z1, z2);
		y += p0 / 2;
	}

	// step
	ofd_geometry(1, 1, x2, x2, y2, y2 - ystep, z2, z2 + dh);
}