/*
Multiloop_modified.c / OpenFDTD / ch.5
*/

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

int main(void)
{
	const int N = 3;   // = 3, 5, 7
	const double fstart = 2.2e9;
	const double fstop  = 3.8e9;
	const int fdiv = 160;
	const double L3 = 100e-3;
	const double C1 = 1.19 * L3;
	const double Cn = 0.87 * L3;
	const double h = 0.1 * L3;
	const double Lv = 0.06 * L3;
	const double Lh = 0.19 * L3;
	const double feed = (C1 + Cn) / 2 / 8;
	const double GP = 0.6 * L3;
	const double d = 0.01 * L3;

	char title[BUFSIZ], fout[BUFSIZ];
	sprintf(title, "modified Multiloop antenna N=%d", N);
	sprintf(fout, "Multiloop_modified_N%d.ofd", N);

	// initialize

	ofd_init();

	// title

	ofd_title(title);

	// mesh

	int div1 = NINT((C1 - Cn) / 8, d);
	if (div1 < N - 1) div1 = N - 1;
	if (div1 % 2 == 1) div1++;

	int div2 = NINT(Cn / 4, d);
	if (div2 % 2 == 1) div2++;

	const double x1 = -GP / 2;
	const double x2 = -C1 / 8;
	const double x3 = -Cn / 8;
	const double x4 = +Cn / 8;
	const double x5 = +C1 / 8;
	const double x6 = +GP / 2;
	ofd_xsection(6, x1, x2, x3, x4, x5, x6);
	ofd_xdivision(5, NINT(x2 - x1, d), div1, div2, div1, NINT(x6 - x5, d));

	const double y1 = -GP / 2;
	const double y2 = -C1 / 8;
	const double y3 = -Cn / 8;
	const double y4 = +Cn / 8;
	const double y5 = +C1 / 8;
	const double y6 = +GP / 2;
	ofd_ysection(6, y1, y2, y3, y4, y5, y6);
	ofd_ydivision(5, NINT(y2 - y1, d), div1, div2, div1, NINT(y6 - y5, d));

	const double z1 = -2 * d;
	const double z2 = 2 * h;
	ofd_zsection(2, z1, z2);
	ofd_zdivision(1, NINT(z2 - z1, d));

	// geometry

	// ground
	ofd_geometry(1, 1, x1, x6, y1, y6, z1, 0);

	// loops
	for (int i = 1; i <= N; i++) {
		const double a = (N > 1) ? ((N - i) * C1 + (i - 1) * Cn) / (N - 1) / 8 : feed;
		ofd_geometry(1, 1, -a, +a, -a, -a, h, h);
		ofd_geometry(1, 1, +a, +a, -a, +a, h, h);
		ofd_geometry(1, 1, -a, +a, +a, +a, h, h);
		ofd_geometry(1, 1, -a, -a, -a, +a, h, h);
	}

	// corners

	ofd_geometry(1, 1, x2, x3, y2, y3, h, h);
	ofd_geometry(1, 1, x2, x3, y4, y5, h, h);
	ofd_geometry(1, 1, x4, x5, y2, y3, h, h);
	ofd_geometry(1, 1, x4, x5, y4, y5, h, h);

/*
	// x2 : 変わらない
	ofd_geometry(1, 1, x2, (x2 + x3) / 2, y2, (y2 + y3) / 2, h, h);
	ofd_geometry(1, 1, x3, (x2 + x3) / 2, y3, (y2 + y3) / 2, h, h);
	ofd_geometry(1, 1, x4, (x4 + x5) / 2, y3, (y2 + y3) / 2, h, h);
	ofd_geometry(1, 1, x5, (x4 + x5) / 2, y2, (y2 + y3) / 2, h, h);
	ofd_geometry(1, 1, x3, (x2 + x3) / 2, y4, (y4 + y5) / 2, h, h);
	ofd_geometry(1, 1, x2, (x2 + x3) / 2, y5, (y4 + y5) / 2, h, h);
	ofd_geometry(1, 1, x4, (x4 + x5) / 2, y4, (y4 + y5) / 2, h, h);
	ofd_geometry(1, 1, x5, (x4 + x5) / 2, y5, (y4 + y5) / 2, h, h);
*/
/*
	// x2 : 変わらない(収束遅い)
	ofd_geometry(1, 1, x2, x3, 0, 0, h, h);
	ofd_geometry(1, 1, x4, x5, 0, 0, h, h);
	ofd_geometry(1, 1, 0, 0, y2, y3, h, h);
	ofd_geometry(1, 1, 0, 0, y4, y5, h, h);
*/
	// InvL
	ofd_geometry(1, 1, feed, feed, 0, 0, 0, Lv);
	ofd_geometry(1, 1, feed, feed, 0, feed, Lv, Lv);
	ofd_geometry(1, 1, feed, 2 * feed - Lh, feed, feed, Lv, Lv);

	// feed

	ofd_feed('Z', feed, 0, 0.1 * d, 1, 0, 50);
	ofd_rfeed(10);

	// frequency

	ofd_frequency1(fstart, fstop, fdiv);

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

	// ABC

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

	// solver

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

	// plot iteration

	ofd_plotiter(1);

	// frequency char.

	ofd_plotref(2, -20, 0, 4);
	ofd_freqdiv(8);

	// output

	ofd_outdata(fout);

	return 0;
}