/*
C_FanSA.c / OpenFDTD / ch.9
*/

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

static double radius(int shape, double y, double xp, double yp, double yc, double r, double ys)
{
	double rad;

	if (y < yp) {
		if (shape == 0) {
			rad = xp * (y + ys) / (yp + ys);  // straight
		}
		else {
			rad = xp * (1 - (y - yp) * (y - yp) / (yp + ys) / (yp + ys));  // parabola
		}
	}
	else {
		//rad = sqrt((r * r) - (y - yc) * (y - yc));  // circle
		rad = xp * sqrt(1 - ((y - yp) * (y - yp)) / ((yc + r - yp) * (yc + r - yp)));  // ellipse
	}

	return rad;
}

int main(void)
{
	const int substrate = 1;   // 0/1:基板なし/あり
	const double fstart =  2.0e9;
	const double fstop  = 11.0e9;
	const int fdiv = 90;
	const double r = 5.44e-3;
	const double yp = 1.5 * r;
	const double xp = sqrt(3) / 2 * r;
	const double GPx = 30e-3;
	const double GPy = 30e-3;
	const double Er = 20;
	const double B = 0.5e-3;
	const double margin = 10e-3;
	const double d = 0.25e-3;           // cell size
	const double dfeed = d;

	const double S = 2 * r;

	const char title[][BUFSIZ] = {"Fan-shaped antenna", "Fan-shaped antenna (substrate)"};
	const char fout[][BUFSIZ] = {"C_FanSA.ofd", "C_FanSA_subst.ofd"};

	// initialize

	ofd_init();

	// title

	ofd_title(title[substrate]);

	// mesh

	const double x0 = - (GPx - (S / 2)) - margin;
	const double x1 = x0 + margin;
	const double x2 = - S / 2;
	const double x3 = 0;
	const double x4 = - x2;
	const double x5 = x4 + margin;
	ofd_xsection(6, x0, x1, x2, x3, x4, x5);
	ofd_xdivision(5, NINT(x1 - x0, d), NINT(x2 - x1, d), NINT(x3 - x2, d), NINT(x4 - x3, d), NINT(x5 - x4, d));

	const double y0 = - dfeed - GPy - margin;
	const double y1 = y0 + margin;
	const double y2 = - dfeed;
	const double y3 = 0;
	const double y4 = y3 + d * NINT(S, d);
	const double y5 = y4 + margin;
	ofd_ysection(6, y0, y1, y2, y3, y4, y5);
	ofd_ydivision(5, NINT(y1 - y0, d), NINT(y2 - y1, d), NINT(y3 - y2, d), NINT(y4 - y3, d), NINT(y5 - y4, d));

	const double z0 = - margin;
	const double z1 = + margin;
	ofd_zsection(2, z0, z1);
	ofd_zdivision(1, NINT(z1 - z0, d));

	// material

	ofd_material(Er, 0, 1, 0, "基板");

	// geometry

	// ground
	ofd_geometry(1, 1, x1, x4, y1, y2, 0, 0);

	// substrate
	if (substrate) {
		ofd_geometry(2, 1, x2, x4, y3, y4, -B, +B);
	}

	// fan
	const int ny = NINT(y4 - y3, d);
	for (int iy = 0; iy < ny; iy++) {
		const double ya = y3 + (iy + 0) * d;
		const double yb = y3 + (iy + 1) * d;
		const double rad = radius(1, (ya + yb) / 2, xp, yp, r, r, dfeed);
		ofd_geometry(1, 1, -rad, +rad, ya, yb, 0, 0);
	}

	// feed
	ofd_geometry(1, 1, x3, x3, y2, y3, 0, 0);

	// feed

	ofd_feed('Y', x3, y2 + 0.1 * d, 0, 1, 0, 50);
	//ofd_rfeed(10);

	// frequency

	ofd_frequency1(fstart, fstop, fdiv);

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

	// ABC

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

	// solver

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

	// plot iteration

	ofd_plotiter(1);

	// frequency char.

	//ofd_plotsmith();
	ofd_plotzin(1, 0, 0, 0);
	//ofd_plotyin(1, 0, 0, 0);
	ofd_plotref(2, -25, 0, 5);
	ofd_freqdiv(9);

	// output

	ofd_outdata(fout[substrate]);

	return 0;
}