/*
PerioAE.c / OpenFDTD / ch.13
*/

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

int main(void)
{
	const double fstart = 7e9;
	const double fstop  = 9e9;
	const int fdiv = 50;
	const int M = 9;
	const int N = 5;
	const double Er1 = 2.6;
	const double Er2 = 2.2;
	const double B1 = 1e-3;
	const double B2 = 1e-3;
	const double d12 = 18e-3;
	const double lpatch = 10.5e-3;
	const double xf = 2.5e-3;
	const double SGP = 135e-3;
	const double sout = 13.5e-3;
	const double ds = 1.5e-3;
	const double w6 = 2.75e-3;
	const double dw = 0.50e-3;
	const double d = 0.25e-3;     // セルサイズ

	// initialize

	ofd_init();

	// title

	ofd_title("Patch Antenna with a Nonuniform Loop Plate");

	// mesh

	const double x0 = - SGP / 2;
	const double x1 = + SGP / 2;
	ofd_xsection(2, x0, x1);
	ofd_xdivision(1, NINT(x1 - x0, d));

	const double y0 = - SGP / 2;
	const double y1 = + SGP / 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 + B1;
	const double z3 = z2 + d12;
	const double z4 = z3 + B2;
	const double z5 = z4 + 15e-3;
	ofd_zsection(6, z0, z1, z2, z3, z4, z5);
	ofd_zdivision(5, NINT(z1 - z0, d), NINT(z2 - z1, d), NINT(z3 - z2, d), NINT(z4 - z3, d), NINT(z5 - z4, d));

	// material

	ofd_material(Er1, 0.0, 1.0, 0.0, "基板1");  // = 2
	ofd_material(Er2, 0.0, 1.0, 0.0, "基板2");  // = 3

	// geometry

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

	// 基板2
	ofd_geometry(3, 1, x0, x1, y0, y1, z3, z4);

	// グラウンド板
	ofd_geometry(1, 1, x0, x1, y0, y1, z0, z1);

	// 給電線
	ofd_geometry(1, 1, xf, xf, 0, 0, z1, z2);

	// パッチアンテナ
	ofd_geometry(1, 1, - lpatch / 2, + lpatch / 2, - lpatch / 2, + lpatch / 2, z2, z2);

	// feed
	ofd_feed('Z', xf, 0, z1 + 0.1 * d, 1, 0, 50);
	ofd_rfeed(10);  // 重要

	// 周期板状ループ
	const double eps = 0.1 * d;
	const double s2 = SGP;
	const double pitch = sout + ds;
	double xc = - s2 / 2 + pitch / 2;
	double yc = xc;
	for (int ix = 0; ix < M; ix++) {
	for (int iy = 0; iy < 2 * N - 1; iy++) {
		const double wm = w6 + (ix - 5) * dw;
		const double x = xc + (ix * pitch);
		const double y = yc + (iy * pitch);
		const double sm = sout - (2 * wm);
		if (wm > -eps) {
			// ループ外周:PEC
			ofd_geometry(1, 1, x - sout / 2, x + sout / 2, y - sout / 2, y + sout / 2, z3, z3);
			// ループ内周:誘電体
			ofd_geometry(3, 1, x - sm / 2 + eps, x + sm / 2 - eps, y - sm / 2 + eps, y + sm / 2 - eps, z3, z3);
		}
	}
	}

	// ABC

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

	// frequency

	ofd_frequency1(fstart, fstop, fdiv);

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

	// solver

	ofd_solver(20000, 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);
	ofd_freqdiv(10);

	// far-1d

	ofd_plotfar1d('Y', 360, 0);
	ofd_plotfar1d('X', 360, 0);

	ofd_far1dcomponent(1, 0, 0);
	ofd_far1dscale(-15, +20, 7);

	// window

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

	// output

	ofd_outdata("PerioAE.ofd");

	return 0;
}