/*
Metaspiral.c / OpenFDTD / ch.22
*/

#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 M = 6;
	const double w = 2.0e-3;
	const double B = 1.6e-3;
	const double er = 2.6;
	const double L1 = 10e-3;
	const double p = L1;
	//const double p0 = 4e-3;
	const double dg = 1e-3;
	//const double rvia = 0.5e-3;
	//const double hvia = 0.6e-3;
	const double Cz = 0.53e-12;
	const double Ly = 3.83e-9;
	const double RB = 80;
	const double margin = 20e-3;
	const double GPx = (M + 1) * L1 + (2 * margin);
	const double GPy = GPx;
	const double d = 1e-3;   // 基本セルサイズ
	char title[BUFSIZ], fout[BUFSIZ];

	sprintf(title, "Metaspiral Antenna M=%d d=%gmm", M, d * 1e3);
	sprintf(fout, "Metaspiral_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 = +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 = 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 dh = B / NINT(B, d);  // 基板中のZ方向セルサイズ
	const double dz = (d + dh) / 2;  // 基板上下の平均
	const double er3 = (2 * Cz) / eps0 * dg / (d * dz);
	ofd_material(er3, 0, 1, 0, "2Cz");

	// geometry

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

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

	// metaspiral
	for (int arm = -1; arm <= 1; arm += 2) {
		double g0 = 0.5 * L1;
		double g1 = 0;
		double g2 = L1;
		for (int m = 0; m < M; m++) {
			if (m % 2 == 0) {
				metaline_y(arm * g0, arm * g1, arm * g2, z1, z2, m + 1, dg, w, d, 3, Ly);
			}
			else {
				metaline_x(arm * g0, arm * g1, arm * g2, z1, z2, m + 1, dg, w, d, 3, Ly);
			}
			g1 = g0;
			g0 = g2;
			g2 = g1 + (m + 2) * L1 * ((m % 4 == 2 || m % 4 == 3) ? +1 : -1);
		}
	}

	// feed + load
	const double xfeed = L1/ 2;
	const double xload = -(M + 1) / 2.0 * p;
	const double yload = +(M + 0) / 2.0 * p;
	for (int arm = -1; arm <= +1; arm += 2) {
		// feed
		ofd_geometry(1, 1, arm * xfeed, arm * xfeed, 0, 0, z1, z2);
		ofd_feed('Z', arm * xfeed, 0, z1 + 0.1 * d, arm, 0, RB);
		// load
		ofd_geometry(1, 1, arm * xload, arm * xload, arm * yload, arm * yload, z1, z2);
		ofd_load('Z', arm * xload, arm * yload, z1 + 0.1 * d, 'R', RB);
	}
	ofd_rfeed(10);

	// ABC

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

	// frequency

	ofd_frequency1(2e9, 4e9, 200);

	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, +15, 7);

	// window

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

	// output

	ofd_outdata(fout);

	return 0;
}

/*
X方向
y0:中心のY座標
x1:始点のX座標
x2:終点のX座標
np:周期の数
*/
static void metaline_x(double y0, double x1, double x2, double z1, double z2,
	int np, double dg, double w, double d, int m, double Ly)
{
	dg = (x2 > x1) ? +dg : -dg;
	const double p = (x2 - x1) / np;
	const double p0 = (p - 2 * dg) / 2;

	//ofd_geometry(1, 1, x1, x2, y0, y0, z2, z2);
	const double y1 = y0 - w / 2;
	const double y2 = y0 + w / 2;

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

		// C
		ofd_geometry(m, 1, x, x + dg, y0, y0, z2, z2);
		//ofd_geometry(1, 1, x, x + dg, y0, y0, z2, z2);
		//ofd_load(dir, x + eps, y0, z2, 'C', 2 * Cz);
		x += dg;

		// patch + via + L
		ofd_geometry(1, 1, x, x + p0, y1, y2, z2, z2);
		ofd_geometry(1, 1, x + p0 / 2, x + p0 / 2, y0, y0, z1, z2);
		ofd_load('Z', x + p0 / 2, y0, z1 + 0.1 * d, 'L', Ly);
		x += p0;

		// C
		ofd_geometry(m, 1, x, x + dg, y0, y0, z2, z2);
		//ofd_geometry(1, 1, x, x + dg, y0, y0, z2, z2);
		//ofd_load(dir, x + eps, y0, z2, 'C', 2 * Cz);
		x += dg;

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

/*
Y方向
x0:中心のX座標
y1:始点のY座標
y2:終点のY座標
np:周期の数
*/
static void metaline_y(double x0, double y1, double y2, double z1, double z2,
	int np, double dg, double w, double d, int m, double Ly)
{
	dg = (y2 > y1) ? +dg : -dg;
	const double p = (y2 - y1) / np;
	const double p0 = (p - 2 * dg) / 2;

	//ofd_geometry(1, 1, x0, x0, y1, y2, z2, z2);
	const double x1 = x0 - w / 2;
	const double x2 = x0 + w / 2;

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

		// C
		ofd_geometry(m, 1, x0, x0, y, y + dg, z2, z2);
		//ofd_geometry(1, 1, x0, x0, y, y + dg, z2, z2);
		//ofd_load(dir, x0, y + eps, z2, 'C', 2 * Cz);
		y += dg;

		// patch + via + L
		ofd_geometry(1, 1, x1, x2, y, y + p0, z2, z2);
		ofd_geometry(1, 1, x0, x0, y + p0 / 2, y + p0 / 2, z1, z2);
		ofd_load('Z', x0, y + p0 / 2, z1 + 0.1 * d, 'L', Ly);
		y += p0;

		// C
		ofd_geometry(m, 1, x0, x0, y, y + dg, z2, z2);
		//ofd_geometry(1, 1, x0, x0, y, y + dg, z2, z2);
		//ofd_load(dir, x0, y + eps, z2, 'C', 2 * Cz);
		y += dg;

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