/*
InvF_EBG.c / OpenFDTD / ch.11
*/

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

int main(void)
{
	const double fstart = 5.5e9;
	const double fstop  = 6.5e9;
	const int fdiv = 50;
	const double L6 = 50e-3;
	const double W = 0.03 * L6;
	const double LV = 0.04 * L6;
	const double B = 0.04 * L6;
	const double Er = 2.2;
	const double DBS = 0.02 * L6;
	const double dfeed = 0.392 * L6;
	const double Spatch = 0.246 * L6;
	const double delta = 0.03 * L6;
	const double margin = 15e-3;
	const double d = 0.5e-3;  // 基本セルサイズ

	// parameters
	const int N = 6;  // even
	const double LH = 0.07 * L6;

	char title[BUFSIZ], fout[BUFSIZ];
	sprintf(title, "InvF-EBG N=%d LH=%.2fL6", N, LH / L6);
	strcpy(fout, "InvF_EBG.ofd");

	const double SEBG = (N * Spatch + (N - 1) * delta);

	// initialize

	ofd_init();

	// title

	ofd_title(title);

	// mesh

	int divpatch = NINT(Spatch, d);
	if (divpatch % 2 == 1) divpatch++;  // even
	int divspace = NINT(delta, d);
	if (divspace % 2 == 1) divspace++;  // even

	double p = - SEBG / 2 - margin;
	ofd_xsection1(p); ofd_ysection1(p);
	// margin
	p += margin;
	ofd_xsection1(p); ofd_ysection1(p);
	ofd_xdivision1(NINT(margin, d)); ofd_ydivision1(NINT(margin, d));
	// EBG
	for (int i = 0; i < N; i++) {
		// patch
		p += Spatch;
		ofd_xsection1(p); ofd_ysection1(p);
		ofd_xdivision1(divpatch); ofd_ydivision1(divpatch);
		// gap
		if (i < N - 1) {
			p += delta;
			ofd_xsection1(p); ofd_ysection1(p);
			ofd_xdivision1(divspace); ofd_ydivision1(divspace);
		}
	}
	// margin
	p += margin;
	ofd_xsection1(p); ofd_ysection1(p);
	ofd_xdivision1(NINT(margin, d)); ofd_ydivision1(NINT(margin, d));

	const double z0 = -2 * d - margin;
	const double z1 = -2 * d;
	const double z2 = 0;
	const double z3 = z2 + B;
	const double z4 = z3 + margin;
	ofd_zsection(5, z0, z1, z2, z3, z4);
	ofd_zdivision(4, NINT(z1 - z0, d), NINT(z2 - z1, d), NINT(z3 - z2, d), NINT(z4 - z3, d));

	// material

	ofd_material(Er, 0, 1, 0, "基板");  // = 2

	// geometry

	// substarate + ground
	const double hgp = SEBG / 2;
	ofd_geometry(2, 1, -hgp, +hgp, -hgp, +hgp, z2, z3);
	ofd_geometry(1, 1, -hgp, +hgp, -hgp, +hgp, z1, z2);

	// EBG
	for (int i = 0; i < N; i++) {
	for (int j = 0; j < N; j++) {
		const double x0 = - hgp + Spatch / 2 + i * (Spatch + delta);
		const double y0 = - hgp + Spatch / 2 + j * (Spatch + delta);
		// patch
		ofd_geometry(1, 1, x0 - Spatch / 2, x0 + Spatch / 2, y0 - Spatch / 2, y0 + Spatch / 2, z3, z3);
		// via
		const double rad = 1 * d;  // 要調整
		ofd_geometry(1, 1, x0 - rad, x0 + rad, y0 - rad, y0 + rad, z2, z3);
	}
	}

	// InvF
	double p0 = delta / 2;
	double p1 = dfeed / 2 - W / 2 - LH;
	double p2 = p1 + LH;
	double p3 = p2 + W / 2;
	double p4 = p3 + W / 2;
	double p5 = p4 + DBS;
	double p6 = p5 + W;
	const double dp = Spatch / divpatch;
	p1 = p0 + dp * NINT(p1 - p0, dp);
	p2 = p0 + dp * NINT(p2 - p0, dp);
	p3 = p0 + dp * NINT(p3 - p0, dp);
	p4 = p0 + dp * NINT(p4 - p0, dp);
	p5 = p0 + dp * NINT(p5 - p0, dp);
	p6 = p0 + dp * NINT(p6 - p0, dp);

	// feed @ -X
	ofd_geometry(1, 1, -p3, -p3, 0, 0, 0, d);

	// InvF : -X and +X
	ofd_geometry(1, 1, -p2, -p4, 0, 0, z2 + d,  z2 + LV    );  // 給電線(ギャップあり)
	ofd_geometry(1, 1, -p5, -p6, 0, 0, z2,      z2 + LV    );  // 短絡線
	ofd_geometry(1, 1, -p1, -p6, 0, 0, z2 + LV, z2 + LV + W);  // 水平放射線
	ofd_geometry(1, 1, +p2, +p4, 0, 0, z2 + d,  z2 + LV    );
	ofd_geometry(1, 1, +p5, +p6, 0, 0, z2,      z2 + LV    );
	ofd_geometry(1, 1, +p1, +p6, 0, 0, z2 + LV, z2 + LV + W);

	// InvF : -Y and +Y
	ofd_geometry(1, 1, 0, 0, -p2, -p4, z2 + d,  z2 + LV    );
	ofd_geometry(1, 1, 0, 0, -p5, -p6, z2,      z2 + LV    );
	ofd_geometry(1, 1, 0, 0, -p1, -p6, z2 + LV, z2 + LV + W);
	ofd_geometry(1, 1, 0, 0, +p2, +p4, z2 + d,  z2 + LV    );
	ofd_geometry(1, 1, 0, 0, +p5, +p6, z2,      z2 + LV    );
	ofd_geometry(1, 1, 0, 0, +p1, +p6, z2 + LV, z2 + LV + W);

	// feed @ -X
	ofd_feed('Z', -p3, 0, z2 + 0.1 * d, 1, 0, 50);
	ofd_rfeed(10);

	// frequency

	ofd_frequency1(fstart, fstop, fdiv);

	const double f0 = 6e9;
	ofd_frequency2(f0, f0, 0);

	// ABC

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

	// solver

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

	// iteration

	ofd_plotiter(1);

	// frequency char.

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

	// output

	ofd_outdata(fout);

	return 0;
}