transducers.h
#ifndef _TRANSDUCERS_H
#define _TRANSDUCERS_H

#include <Eigen/Core>
#include <complex>
#include <cmath>
#include <iostream>

using namespace std;

using Eigen::VectorXd;
using Eigen::MatrixXd;
using Eigen::MatrixXcd;
using Eigen::VectorXcd;

class TransducerLayout {
	VectorXd transX, transY;
	static const double P0 = 5;
	static const double k = 732.733;
	static const double r = 0.005;
public:
	TransducerLayout() {}
	TransducerLayout(VectorXd tx, VectorXd ty) {
		transX = tx;
		transY = ty;
	}
	
	MatrixXcd getPressure(VectorXd x, VectorXd y, VectorXd z) {
		int numTrans = transX.rows();
		int numPP = x.rows();
		MatrixXcd p (numTrans, numPP);

		for (int j = 0; j < numPP; j++) {
			p.col(j) = getP(x(j),y(j),z(j));
		}

		return p;
	}

	VectorXcd getP(double x, double y, double z) {
		int numTrans = transX.rows();
		VectorXcd P (numTrans);

		for (int i = 0; i < numTrans; i++) {
			double xr = x - transX(i);
			double yr = y - transY(i);
			double d = sqrt(pow(xr,2) + pow(yr,2) + pow(z,2));
			double rho = sqrt(pow(xr,2) + pow(yr,2));
			if (rho != 0)
				P(i) = 2.0*complexExp(k*d)*P0*jn(1,k*r*rho/d)/(k*r*rho);
			else
				P(i) = P0*complexExp(k*d)/d;
		}

		return P;
	}

	VectorXcd getPx(double x, double y, double z) {
		int numTrans = transX.rows();
		VectorXcd Px (numTrans);

		for (int i = 0; i < numTrans; i++) {
			double xr = x - transX(i);
			double yr = y - transY(i);
			double d = sqrt(pow(xr,2) + pow(yr,2) + pow(z,2));
			double rho = sqrt(pow(xr,2) + pow(yr,2));
			if (rho != 0)
				Px(i) = (complexExp(d*k)*P0*xr*(k*r*rho*(-pow(d,2) + pow(rho,2))*jn(0, (k*r*rho)/d) + 2.0*pow(d,2)*(d - 1i*k*pow(rho,2))*jn(1, (k*r*rho)/d) + k*r*rho*(pow(d,2) - pow(rho,2))*jn(2, (k*r*rho)/d)))/(2.0*pow(d,3)*k*r*pow(rho,3));
			else
				Px(i) = P0*complexExp(k*d)/d;
		}

		return Px;
	}

	VectorXcd getPy(double x, double y, double z) {
		int numTrans = transX.rows();
		VectorXcd Py (numTrans);

		for (int i = 0; i < numTrans; i++) {
			double xr = x - transX(i);
			double yr = y - transY(i);
			double d = sqrt(pow(xr,2) + pow(yr,2) + pow(z,2));
			double rho = sqrt(pow(xr,2) + pow(yr,2));
			if (rho != 0)
				Py(i) = -((complexExp(d*k)*P0*yr*(k*r*rho*(-pow(d,2) + pow(rho,2))*jn(0, (k*r*rho)/d) + 2.0*pow(d,2)*(d - 1i*k*pow(rho,2))*jn(1, (k*r*rho)/d) + k*r*rho*(pow(d,2) - pow(rho,2))*jn(2, (k*r*rho)/d)))/(2.0*pow(d,3)*k*r*pow(rho,3)));
			else
				Py(i) = P0*complexExp(k*d)/d;
		}

		return Py;
	}

	VectorXcd getPz(double x, double y, double z) {
		int numTrans = transX.rows();
		VectorXcd Pz (numTrans);

		for (int i = 0; i < numTrans; i++) {
			double xr = x - transX(i);
			double yr = y - transY(i);
			double d = sqrt(pow(xr,2) + pow(yr,2) + pow(z,2));
			double rho = sqrt(pow(xr,2) + pow(yr,2));
			if (rho != 0)
				Pz(i) = (complexExp(d*k)*P0*z*((-r)*rho*jn(0, (k*r*rho)/d) + 2.0*pow(d,2)*1i*jn(1, (k*r*rho)/d) + r*rho*jn(2, (k*r*rho)/d)))/(2.0*pow(d,3)*r*rho);
			else
				Pz(i) = P0*complexExp(k*d)/d;
		}

		return Pz;
	}

	VectorXcd getPxx(double x, double y, double z) {
		int numTrans = transX.rows();
		VectorXcd Pxx (numTrans);

		for (int i = 0; i < numTrans; i++) {
			double xr = x - transX(i);
			double yr = y - transY(i);
			double d = sqrt(pow(xr,2) + pow(yr,2) + pow(z,2));
			double rho = sqrt(pow(xr,2) + pow(yr,2));
			if (rho != 0)
				Pxx(i) = (1.0/(2.0*r*pow(rho,5)))*(complexExp(d*k)*P0*((r*rho*(pow(d,2) - pow(rho,2))*(-3.0*pow(rho,2)*pow(xr,2) + 2.0*d*1i*k*pow(rho,2)*pow(xr,2) + pow(d,2)*(pow(rho,2) - 3.0*pow(xr,2)))*jn(0, (k*r*rho)/d))/pow(d,5) + (1.0/2)*(((k*pow(rho,2)*(6.0*pow(d,2)*pow(r,2)*pow(rho,2) - 3.0*pow(r,2)*pow(rho,4) + pow(d,4)*(-3.0*pow(r,2) - 4.0*pow(rho,2)))*pow(xr,2))/pow(d,6) - (4.0*(pow(rho,2) - 3.0*pow(xr,2)))/k + (4.0*1i*pow(rho,2)*((-pow(rho,2))*pow(xr,2) + pow(d,2)*(pow(rho,2) - 2.0*pow(xr,2))))/pow(d,3))*jn(1, (k*r*rho)/d) + (r*rho*(-pow(d,2) + pow(rho,2))*(2.0*d*(-3.0*pow(rho,2)*pow(xr,2) + 2.0*d*1i*k*pow(rho,2)*pow(xr,2) +pow(d,2)*(pow(rho,2) - 3.0*pow(xr,2)))*jn(2, (k*r*rho)/d) + k*r*rho*(-pow(d,2) + pow(rho,2))*pow(xr,2)*jn(3, (k*r*rho)/d)))/pow(d,6))));
			else
				Pxx(i) = P0*complexExp(k*d)/d;
		}

		return Pxx;
	}

	VectorXcd getPxy(double x, double y, double z) {
		int numTrans = transX.rows();
		VectorXcd Pxy (numTrans);

		for (int i = 0; i < numTrans; i++) {
			double xr = x - transX(i);
			double yr = y - transY(i);
			double d = sqrt(pow(xr,2) + pow(yr,2) + pow(z,2));
			double rho = sqrt(pow(xr,2) + pow(yr,2));
			if (rho != 0)
				Pxy(i) = -((1.0/(4.0*pow(d,6)*k*r*pow(rho,5)))*(complexExp(d*k)*P0*xr*yr*(2.0*d*k*r*rho*(-3.0*pow(d,4) + 2.0*pow(d,3)*1i*k*pow(rho,2) + 3.0*pow(rho,4) - 2.0*d*1i*k*pow(rho,4))*jn(0, (k*r*rho)/d) + (12.0*pow(d,6) - 8.0*pow(d,5)*1i*k*pow(rho,2) - 4.0*pow(d,3)*1i*k*pow(rho,4) + 6.0*pow(d,2)*pow(k,2)*pow(r,2)*pow(rho,4) - 3.0*pow(k,2)*pow(r,2)*pow(rho,6) + pow(d,4)*pow(k,2)*pow(rho,2)*(-3.0*pow(r,2) - 4.0*pow(rho,2)))*jn(1, (k*r*rho)/d) + k*r*rho*(pow(d,2) - pow(rho,2))*(2.0*d*(3.0*pow(d,2) + 3.0*pow(rho,2) - 2.0*d*1i*k*pow(rho,2))*jn(2, (k*r*rho)/d) + k*r*rho*(pow(d,2) - pow(rho,2))*jn(3, (k*r*rho)/d)))));
			else
				Pxy(i) = P0*complexExp(k*d)/d;
		}

		return Pxy;
	}

	VectorXcd getPxz(double x, double y, double z) {
		int numTrans = transX.rows();
		VectorXcd Pxz (numTrans);

		for (int i = 0; i < numTrans; i++) {
			double xr = x - transX(i);
			double yr = y - transY(i);
			double d = sqrt(pow(xr,2) + pow(yr,2) + pow(z,2));
			double rho = sqrt(pow(xr,2) + pow(yr,2));
			if (rho != 0)
				Pxz(i) = (1.0/(4.0*pow(d,6)*r*pow(rho,3)))*(complexExp(d*k)*P0*xr*z*(-2.0*d*r*rho*(pow(d,3)*1i*k + 3.0*pow(rho,2) - 2.0*d*1i*k*pow(rho,2))*jn(0, (k*r*rho)/d) + (4.0*pow(d,5)*1i + 4.0*pow(d,3)*1i*pow(rho,2) + 4.0*pow(d,4)*k*pow(rho,2) - 3.0*pow(d,2)*k*pow(r,2)*pow(rho,2) + 3.0*k*pow(r,2)*pow(rho,4))*jn(1, (k*r*rho)/d) + r*rho*(2.0*d*(pow(d,3)*1i*k + 3.0*pow(rho,2) - 2.0*d*1i*k*pow(rho,2))*jn(2, (k*r*rho)/d) + k*r*rho*(pow(d,2) - pow(rho,2))*jn(3, (k*r*rho)/d))));
			else
				Pxz(i) = P0*complexExp(k*d)/d;
		}

		return Pxz;
	}

	VectorXcd getPyy(double x, double y, double z) {
		int numTrans = transX.rows();
		VectorXcd Pyy (numTrans);

		for (int i = 0; i < numTrans; i++) {
			double xr = x - transX(i);
			double yr = y - transY(i);
			double d = sqrt(pow(xr,2) + pow(yr,2) + pow(z,2));
			double rho = sqrt(pow(xr,2) + pow(yr,2));
			if (rho != 0)
				Pyy(i) = -((1.0/(4.0*pow(d,6)*k*r*pow(rho,5)))*(complexExp(d*k)*P0*(-2.0*d*k*r*rho*(pow(d,2) - pow(rho,2))*(-3.0*pow(rho,2)*pow(yr,2) + 2.0*d*1i*k*pow(rho,2)*pow(yr,2) + pow(d,2)*(pow(rho,2) - 3.0*pow(yr,2)))*jn(0, (k*r*rho)/d) + (4.0*pow(d,3)*1i*k*pow(rho,4)*pow(yr,2) - 6.0*pow(d,2)*pow(k,2)*pow(r,2)*pow(rho,4)*pow(yr,2) + 3.0*pow(k,2)*pow(r,2)*pow(rho,6)*pow(yr,2) + pow(d,4)*pow(k,2)*pow(rho,2)*(3.0*pow(r,2) + 4.0*pow(rho,2))*pow(yr,2) + 4.0*pow(d,6)*(pow(rho,2) - 3.0*pow(yr,2)) -4.0*pow(d,5)*1i*k*pow(rho,2)*(pow(rho,2) - 2.0*pow(yr,2)))*jn(1, (k*r*rho)/d) + k*r*rho*(pow(d,2) - pow(rho,2))*(2.0*d*(-3.0*pow(rho,2)*pow(yr,2) + 2.0*d*1i*k*pow(rho,2)*pow(yr,2) + pow(d,2)*(pow(rho,2) - 3.0*pow(yr,2)))*jn(2, (k*r*rho)/d) + k*r*rho*(-pow(d,2) + pow(rho,2))*pow(yr,2)*jn(3, (k*r*rho)/d)))));
			else
				Pyy(i) = P0*complexExp(k*d)/d;
		}

		return Pyy;
	}

	VectorXcd getPyz(double x, double y, double z) {
		int numTrans = transX.rows();
		VectorXcd Pyz (numTrans);

		for (int i = 0; i < numTrans; i++) {
			double xr = x - transX(i);
			double yr = y - transY(i);
			double d = sqrt(pow(xr,2) + pow(yr,2) + pow(z,2));
			double rho = sqrt(pow(xr,2) + pow(yr,2));
			if (rho != 0)
				Pyz(i) = -((1.0/(4.0*pow(d,6)*r*pow(rho,3)))*(complexExp(d*k)*P0*yr*z*(-2.0*d*r*rho*(pow(d,3)*1i*k + 3.0*pow(rho,2) - 2.0*d*1i*k*pow(rho,2))*jn(0, (k*r*rho)/d) + (4.0*pow(d,5)*1i + 4.0*pow(d,3)*1i*pow(rho,2) + 4.0*pow(d,4)*k*pow(rho,2) - 3.0*pow(d,2)*k*pow(r,2)*pow(rho,2) + 3.0*k*pow(r,2)*pow(rho,4))*jn(1, (k*r*rho)/d) + r*rho*(2.0*d*(pow(d,3)*1i*k + 3.0*pow(rho,2) - 2.0*d*1i*k*pow(rho,2))*jn(2, (k*r*rho)/d) + k*r*rho*(pow(d,2) - pow(rho,2))*jn(3, (k*r*rho)/d)))));
			else
				Pyz(i) = P0*complexExp(k*d)/d;;
		}

		return Pyz;
	}

	__complex__ double complexExp(double x) {
		return cos(x) + 1j*sin(x);
	}

	unsigned int getNumberOfTransducers() {
		return transX.rows();
	}
};

#endif