/*
readout.c

read oth.out
*/

#include "oth.h"

void readout(FILE *fp)
{
	size_t num = 0;

	num += fread(Title,       sizeof(char),      256, fp);
	num += fread(&Nx,         sizeof(int),         1, fp);
	num += fread(&Ny,         sizeof(int),         1, fp);
	num += fread(&Nz,         sizeof(int),         1, fp);
	num += fread(&NE,         sizeof(int64_t),     1, fp);
	num += fread(&NH,         sizeof(int64_t),     1, fp);
	num += fread(&NEXI,       sizeof(int64_t),     1, fp);
	num += fread(&NEXJ,       sizeof(int64_t),     1, fp);
	num += fread(&NEXK,       sizeof(int64_t),     1, fp);
	num += fread(&NEX0,       sizeof(int64_t),     1, fp);
	num += fread(&NEYI,       sizeof(int64_t),     1, fp);
	num += fread(&NEYJ,       sizeof(int64_t),     1, fp);
	num += fread(&NEYK,       sizeof(int64_t),     1, fp);
	num += fread(&NEY0,       sizeof(int64_t),     1, fp);
	num += fread(&NEZI,       sizeof(int64_t),     1, fp);
	num += fread(&NEZJ,       sizeof(int64_t),     1, fp);
	num += fread(&NEZK,       sizeof(int64_t),     1, fp);
	num += fread(&NEZ0,       sizeof(int64_t),     1, fp);
	num += fread(&NHXI,       sizeof(int64_t),     1, fp);
	num += fread(&NHXJ,       sizeof(int64_t),     1, fp);
	num += fread(&NHXK,       sizeof(int64_t),     1, fp);
	num += fread(&NHX0,       sizeof(int64_t),     1, fp);
	num += fread(&NHYI,       sizeof(int64_t),     1, fp);
	num += fread(&NHYJ,       sizeof(int64_t),     1, fp);
	num += fread(&NHYK,       sizeof(int64_t),     1, fp);
	num += fread(&NHY0,       sizeof(int64_t),     1, fp);
	num += fread(&NHZI,       sizeof(int64_t),     1, fp);
	num += fread(&NHZJ,       sizeof(int64_t),     1, fp);
	num += fread(&NHZK,       sizeof(int64_t),     1, fp);
	num += fread(&NHZ0,       sizeof(int64_t),     1, fp);
	num += fread(&NMaterial,  sizeof(int),         1, fp);
	num += fread(&NFeed,      sizeof(int),         1, fp);
	num += fread(&IPlanewave, sizeof(int),         1, fp);
	num += fread(&NPoint,     sizeof(int),         1, fp);
	num += fread(&NFrequency, sizeof(int),         1, fp);
	num += fread(&Planewave,  sizeof(planewave_t), 1, fp);
	num += fread(&Pbc.xpbc,   sizeof(int),         1, fp);
	num += fread(&Pbc.ypbc,   sizeof(int),         1, fp);
	num += fread(&Pbc.zpbc,   sizeof(int),         1, fp);
	num += fread(&NGline,     sizeof(int),         1, fp);

	Xn  = (double *)malloc((Nx + 1) * sizeof(double));
	Yn  = (double *)malloc((Ny + 1) * sizeof(double));
	Zn  = (double *)malloc((Nz + 1) * sizeof(double));
	Xc  = (double *)malloc((Nx + 0) * sizeof(double));
	Yc  = (double *)malloc((Ny + 0) * sizeof(double));
	Zc  = (double *)malloc((Nz + 0) * sizeof(double));
	DXn = (double *)malloc((Nx + 1) * sizeof(double));
	DYn = (double *)malloc((Ny + 1) * sizeof(double));
	DZn = (double *)malloc((Nz + 1) * sizeof(double));
	DXc = (double *)malloc((Nx + 0) * sizeof(double));
	DYc = (double *)malloc((Ny + 0) * sizeof(double));
	DZc = (double *)malloc((Nz + 0) * sizeof(double));
	Material  =     (material_t *)malloc(NMaterial * NFrequency * sizeof(material_t));
	Feed      =         (feed_t *)malloc(NFeed                  * sizeof(feed_t));
	Point     =        (point_t *)malloc(NPoint                 * sizeof(point_t));
	Frequency =         (double *)malloc(NFrequency             * sizeof(double));
	NIter     =            (int *)malloc(NFrequency             * sizeof(int));
	Zin       =    (d_complex_t *)malloc(NFeed * NFrequency     * sizeof(d_complex_t));
	Ref       =         (double *)malloc(NFeed * NFrequency     * sizeof(double));
	Spara     =    (d_complex_t *)malloc(NPoint * NFrequency    * sizeof(d_complex_t));
	PLoss     =         (double *)malloc(NFrequency             * sizeof(double));
	Gline     = (double (*)[2][3])malloc(NGline * 2 * 3         * sizeof(double));

	num += fread(Xn,          sizeof(double),      Nx + 1,                 fp);
	num += fread(Yn,          sizeof(double),      Ny + 1,                 fp);
	num += fread(Zn,          sizeof(double),      Nz + 1,                 fp);
	num += fread(Xc,          sizeof(double),      Nx,                     fp);
	num += fread(Yc,          sizeof(double),      Ny,                     fp);
	num += fread(Zc,          sizeof(double),      Nz,                     fp);
	num += fread(DXn,         sizeof(double),      Nx + 1,                 fp);
	num += fread(DYn,         sizeof(double),      Ny + 1,                 fp);
	num += fread(DZn,         sizeof(double),      Nz + 1,                 fp);
	num += fread(DXc,         sizeof(double),      Nx,                     fp);
	num += fread(DYc,         sizeof(double),      Ny,                     fp);
	num += fread(DZc,         sizeof(double),      Nz,                     fp);
	num += fread(Material,    sizeof(material_t),  NFrequency * NMaterial, fp);
	num += fread(Feed,        sizeof(feed_t),      NFeed,                  fp);
	num += fread(Point,       sizeof(point_t),     NPoint,                 fp);
	num += fread(Frequency,   sizeof(double),      NFrequency,             fp);
	num += fread(NIter,       sizeof(int),         NFrequency,             fp);
	num += fread(Zin,         sizeof(d_complex_t), NFrequency * NFeed,     fp);
	num += fread(Ref,         sizeof(double),      NFrequency * NFeed,     fp);
	num += fread(Spara,       sizeof(d_complex_t), NFrequency * NPoint,    fp);
	num += fread(PLoss,       sizeof(double),      NFrequency,             fp);
	num += fread(Gline,       sizeof(double),      NGline * 2 * 3,         fp);

	Residual = (double **)malloc(NFrequency * sizeof(double *));
	Eaverage = (double **)malloc(NFrequency * sizeof(double *));
	for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
		Residual[ifreq] = (double *)malloc(NIter[ifreq] * sizeof(double));
		Eaverage[ifreq] = (double *)malloc(NIter[ifreq] * sizeof(double));
		num += fread(Residual[ifreq], sizeof(double), NIter[ifreq], fp);
		num += fread(Eaverage[ifreq], sizeof(double), NIter[ifreq], fp);
	}

	mE = (int *)malloc(NE * sizeof(int));
	mH = (int *)malloc(NH * sizeof(int));
	num += fread(mE, sizeof(int), NE, fp);
	num += fread(mH, sizeof(int), NH, fp);

	Ex = (f_complex_t ****)malloc(NFrequency * sizeof(f_complex_t ***));
	Ey = (f_complex_t ****)malloc(NFrequency * sizeof(f_complex_t ***));
	Ez = (f_complex_t ****)malloc(NFrequency * sizeof(f_complex_t ***));
	for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
		Ex[ifreq] = (f_complex_t ***)malloc((Nx + 0) * sizeof(f_complex_t **));
		for (int i = 0; i < Nx + 0; i++) {
			Ex[ifreq][i] = (f_complex_t **)malloc((Ny + 1) * sizeof(f_complex_t *));
			for (int j = 0; j < Ny + 1; j++) {
				Ex[ifreq][i][j] = (f_complex_t *)malloc((Nz + 1) * sizeof(f_complex_t));
				num += fread(Ex[ifreq][i][j], sizeof(f_complex_t), Nz + 1, fp);
			}
		}
		Ey[ifreq] = (f_complex_t ***)malloc((Nx + 1) * sizeof(f_complex_t **));
		for (int i = 0; i < Nx + 1; i++) {
			Ey[ifreq][i] = (f_complex_t **)malloc((Ny + 0) * sizeof(f_complex_t *));
			for (int j = 0; j < Ny + 0; j++) {
				Ey[ifreq][i][j] = (f_complex_t *)malloc((Nz + 1) * sizeof(f_complex_t));
				num += fread(Ey[ifreq][i][j], sizeof(f_complex_t), Nz + 1, fp);
			}
		}
		Ez[ifreq] = (f_complex_t ***)malloc((Nx + 1) * sizeof(f_complex_t **));
		for (int i = 0; i < Nx + 1; i++) {
			Ez[ifreq][i] = (f_complex_t **)malloc((Ny + 1) * sizeof(f_complex_t *));
			for (int j = 0; j < Ny + 1; j++) {
				Ez[ifreq][i][j] = (f_complex_t *)malloc((Nz + 0) * sizeof(f_complex_t));
				num += fread(Ez[ifreq][i][j], sizeof(f_complex_t), Nz + 0, fp);
			}
		}
	}

	size_t num0;
	size_t size = fread(&num0, sizeof(size_t), 1, fp);
	size = size;  // suppress gcc warning
	//printf("%zd %zd\n", num, num0);

	if (num != num0) {
		printf("*** invalid file length : (%zd, %zd)\n", num0, num);
	}
}
