Java Vadose Zone

package br.com.sed.seris.domain.simulation.concentration;

import java.io.Serializable;
import java.math.BigDecimal;
import java.math.MathContext;
import java.math.RoundingMode;
import java.text.DecimalFormat;
import java.text.NumberFormat;
import java.util.ArrayList;
import java.util.List;

import org.apache.commons.math3.special.Erf;

import br.com.sed.seris.domain.Chemical;
import br.com.sed.seris.domain.Soil;
import br.com.sed.seris.domain.simulation.Simulation;

public class VadoseZone implements Serializable {
	
	/**
	 * 
	 */
	private static final long serialVersionUID = 1L;
	
	private static  BigDecimal SQRT_DIG = new BigDecimal(150);
	private static final BigDecimal SQRT_PRE = new BigDecimal(10).pow(SQRT_DIG.intValue());
	
	private static final BigDecimal VADOSE_ZONE = new BigDecimal("10"); //TODO: Entradas de usuário;
	private static final BigDecimal DEPTH_TO_TOP = new BigDecimal("2");
	private static final BigDecimal THICKNESS_OF_SOURCE = new BigDecimal("1");
	private static final BigDecimal LENGTH_OF_SOURCE = new BigDecimal("10");
	private static final BigDecimal WIDTH_OF_SOURCE = new BigDecimal("10");
	
	
	
	public Double getConcentration(Simulation simulation) {
		
		for (int i = 0; i < simulation.getSimulationChemicals().size(); i++) {
			//PASSO 1.0.0:
			BigDecimal si = new BigDecimal(simulation.getSimulationChemicals().get(i).getChemical().getSolubility()).setScale(30, RoundingMode.HALF_EVEN);
			BigDecimal cti =  new BigDecimal("0.000016").setScale(30, RoundingMode.HALF_EVEN);
			BigDecimal mwTPH = new BigDecimal("80").setScale(30, RoundingMode.HALF_EVEN);
			BigDecimal mwi = new BigDecimal(simulation.getSimulationChemicals().get(i).getChemical().getMolecularWeight()).setScale(30, RoundingMode.HALF_EVEN);
			BigDecimal cTPH = new BigDecimal("0.00002").setScale(30, RoundingMode.HALF_EVEN); //(TPH)
			BigDecimal selfi = BigDecimal.ZERO;
			selfi =  (cti.divide(cTPH, MathContext.DECIMAL128)).multiply(mwTPH.divide(mwi, MathContext.DECIMAL128)).multiply(si).setScale(30, RoundingMode.HALF_EVEN);
			System.out.println("Valor de selfi :D");	
		}
		
		
		
		//simulation.getSoil().getSoilBulkDensity()
		
		
		return null;
	}
	
	
	public static void main(String[] args) {
		
		NumberFormat nb = new DecimalFormat("0.##E0");
		Chemical chemical = new Chemical();
		chemical.setSolubility((float) 1.75E+03);
		chemical.setMolecularWeight((float) 78);
		
		// Passo 1.0.0 : cálculo de concentração inicial da fase dissolvida
		BigDecimal si = new BigDecimal(chemical.getSolubility()).setScale(30, RoundingMode.HALF_EVEN);
		BigDecimal cti =  new BigDecimal("0.000016").setScale(30, RoundingMode.HALF_EVEN);
		BigDecimal mwTPH = new BigDecimal("80").setScale(30, RoundingMode.HALF_EVEN);
		BigDecimal mwi = new BigDecimal(chemical.getMolecularWeight()).setScale(30, RoundingMode.HALF_EVEN);
		BigDecimal cTPH = new BigDecimal("0.00002").setScale(30, RoundingMode.HALF_EVEN); //(TPH)
		BigDecimal selfi = BigDecimal.ZERO;
		selfi =  (cti.divide(cTPH, MathContext.DECIMAL128)).multiply(mwTPH.divide(mwi, MathContext.DECIMAL128)).multiply(si).setScale(30, RoundingMode.HALF_EVEN);
		
		// Passo 1.2.1 : Cálculo da porosidade preenchida de água na zona vadosa
		Soil soil = new Soil();
		soil.setHydraulicConductivity(25f);
		soil.setResidualWaterContent(0.21f);
		soil.setTotalPorosity(0.35f);
		BigDecimal q = new BigDecimal("0.054794520547945200000000000000").setScale(30, RoundingMode.HALF_EVEN); //Entrada do usuário 
		BigDecimal Ksuv = new BigDecimal(soil.getHydraulicConductivity()).setScale(30, RoundingMode.HALF_EVEN);
		BigDecimal KsuLens = new BigDecimal(soil.getHydraulicConductivity()).setScale(30, RoundingMode.HALF_EVEN);
		BigDecimal krv = BigDecimal.ZERO;
		BigDecimal Orv = new BigDecimal(soil.getResidualWaterContent().toString()).setScale(30, RoundingMode.HALF_EVEN);
		BigDecimal Otv = new BigDecimal(soil.getTotalPorosity().toString()).setScale(30, RoundingMode.HALF_EVEN);
		BigDecimal Nl = new BigDecimal("2.68").setScale(30, RoundingMode.HALF_EVEN); 
		BigDecimal Nv = new BigDecimal("1.37").setScale(30, RoundingMode.HALF_EVEN); //TODO: feature in preview
		
		if (q.compareTo(KsuLens) >= 0 ) {
			KsuLens = q;
		}
		krv = q.divide(KsuLens).setScale(30, RoundingMode.HALF_EVEN);
		
		BigDecimal Yv = calcularParametroDeDistribuicaoDeTamanhoDeProZonaVadosa(Nv);
		
		BigDecimal cm3 = new BigDecimal("1");
		BigDecimal Owv = calcularPorosidadeDaAgua(krv, Orv, Otv, Yv);
		
		// Passo 1.2.2 : Cálculo da concentração da fase dissolvida (using equilibrium partitioning - 3 phases) (Cw)
		BigDecimal Oav = calcularPorosidadeDoArZonaLens(Otv, Owv);
		BigDecimal KOC =  new BigDecimal("59").setScale(30, RoundingMode.HALF_EVEN); //TODO: declaracao
		BigDecimal gm = new BigDecimal("1000").setScale(30, RoundingMode.HALF_EVEN);
		BigDecimal L = new BigDecimal("1000").setScale(30, RoundingMode.HALF_EVEN);
		BigDecimal mg =  new BigDecimal("1000").setScale(30, RoundingMode.HALF_EVEN);
		BigDecimal pb = new BigDecimal("1.7").setScale(30, RoundingMode.HALF_EVEN);
		BigDecimal koc =  new BigDecimal("59").multiply(cm3.divide(gm)).setScale(30, RoundingMode.HALF_EVEN);
		BigDecimal foc = new BigDecimal("0.005").multiply(gm.divide(gm)).setScale(30, RoundingMode.HALF_EVEN);
		double khPow = Math.pow(10, -1);
		BigDecimal kh = new BigDecimal("2.28").multiply(new BigDecimal(khPow)).multiply(mg.divide(L).divide(mg.divide(L))).setScale(30, RoundingMode.HALF_EVEN);		BigDecimal pbGL = pb.multiply(new BigDecimal(Math.pow(10,3))).setScale(30, RoundingMode.HALF_EVEN);
		BigDecimal cw = calcularConcentracaoInicialFaseDisolvida(cti, Owv, Oav, koc,foc, kh, pbGL);
		
		//PASSO 1.2.3 : Concentração inicial na fonte (Cw0)
		BigDecimal cwO = calcularConcentracaoInicialDaFonte(si, selfi, cw);
		
		// Passo 2 : Cálculo da concentração de vapor Inicial (Cv0)
		calcularConcentracaoDeVaporInicial(kh, cwO);
				
		// PASSO 3 :  CÁLCULO DE BETA (β)
		BigDecimal Krl = calcularPermeabilidadeRelativa(q);
				
		BigDecimal Orl = new BigDecimal("0.05").multiply(cm3.divide(cm3)).setScale(30, RoundingMode.HALF_EVEN); //TODO: Declaracao 
		BigDecimal Otl = new BigDecimal("0.30").multiply(cm3.divide(cm3)).setScale(30, RoundingMode.HALF_EVEN);

		BigDecimal yL = calcularParametroDeDistribuicaoTamanhoDoPoroZonaLens(Nl);
		
		BigDecimal Owl = calcularPorosidadeDaAgua(Krl, Orl, Otl, yL);
		BigDecimal Oal = calcularPorosidadeDoArZonaLens(Otl, Owl);
		
		BigDecimal s = new BigDecimal("1"); //TODO: declaracao
		BigDecimal cm2 = new BigDecimal("1");
		BigDecimal Dw = new BigDecimal("9.8").multiply(new BigDecimal(Math.pow(10, -6)).multiply(cm2.divide(s, MathContext.DECIMAL128)));
		BigDecimal Da = new BigDecimal("8.8").multiply(new BigDecimal(Math.pow(10, -2)).multiply(cm2.divide(s, MathContext.DECIMAL128)));
		
		BigDecimal Deffv = calcularCoeficienteDeDifusaoEfetivaZonaVadosa(Otv, Owv, Oav, kh, Dw, Da);
		BigDecimal Deffl = calcularCoeficienteDeDifusaoEfetivaZonaLens(kh, Otl,	Owl, Oal, Dw, Da);

		BigDecimal L1 = new BigDecimal("0"); //m //TODO:declaracao
		BigDecimal Lv = new BigDecimal("10"); //m
		BigDecimal Lv1 = Lv.subtract(L1);
		BigDecimal Defft = calcularCoeficienteDeDifusaoEfetivaTotal(Deffv, Deffl, L1, Lv, Lv1);
	    
		BigDecimal Lw = new BigDecimal("1"); //TODO: declaracao
		BigDecimal siGcm3 = conveterSiParaCM3(si);
		BigDecimal Psource = new BigDecimal("2"); //TODO: declaração 

		BigDecimal Ld = Psource.add(Lw.divide(new BigDecimal("2")));
		BigDecimal defftCm2D = converterDefftParaCm2PorDia(Defft);
		BigDecimal LwCM = converterLwParaCm(Lw);
		Ld = calcularDistanciaSuperficieDoSoloCentroDaFonte(Lw, Psource);
		
		System.out.println("valor de ld:" + Ld);
		System.out.println("valor de lw:" + Lw);
		System.out.println("valor de Psource:" + Psource);
		
		BigDecimal Bcrp = BigDecimal.ZERO;
		BigDecimal Bcrv = BigDecimal.ZERO;
		BigDecimal Bcr = BigDecimal.ZERO;
		BigDecimal Bsrp = BigDecimal.ZERO;
		BigDecimal Bsrv = BigDecimal.ZERO;
		BigDecimal Bsr = BigDecimal.ZERO;
		
		/**
		 * Verifica se tem fase residual
		 */
		if (cw.compareTo(selfi) >=0 ) {
			Bcrp = calcularQuedaDaFontePorPercolacaoComFaseResidual(mwTPH, mwi,	cTPH, q, pb, siGcm3, LwCM);
			Bcrv = calcularQuedaDaFontePorVolatizacaoComFaseResidual(mwTPH, mwi, cTPH, pb, kh, siGcm3, Ld, defftCm2D, LwCM);
			Bcr = calcularQuedaDaFonteSemFaseResidual(Bcrp, Bcrv);
		} else {
			Bsrp = calcularQuedaDaFontePorPercolacaoSemFaseResidual(nb, q, Owv, Oav, KOC, pb, foc, kh, LwCM);
			Bsrv = calcularQuedaDaFontePorVolatizacaoSemFaseResidual(nb, Owv, Oav, KOC, pb, foc, kh, Ld, defftCm2D, LwCM);
			Bsr = calcularQuedaDaFonteSemFaseResidual(Bsrp, Bsrv);
		}
		
		BigDecimal t = new BigDecimal("365").setScale(30, RoundingMode.HALF_EVEN);
		List<BigDecimal> listBxt = new ArrayList<BigDecimal>();
		BigDecimal xSource = LENGTH_OF_SOURCE.multiply(new BigDecimal("100"));
		BigDecimal ySource = WIDTH_OF_SOURCE.multiply(new BigDecimal("100"));
		
		// Passo 4 : CÁLCULO DE B (x,t)
		BigDecimal v = calcularVelocidadeDePercolacao(q, Owv);
		BigDecimal Xm = calcularDistanciaDaFonteAteLencolFreatico(Lv, Lw,	Psource);
		BigDecimal Xm1 = new BigDecimal("1");
		BigDecimal XmO = Xm.divide(Xm1);
		BigDecimal aLx = calcularDoCoeficienteDeDispersao(Xm);
		BigDecimal Dx = aLx.multiply(v).multiply(new BigDecimal("100"));
		BigDecimal R = calcularFatorDeRetardacao(Owv, KOC, pb, foc);
		BigDecimal W = calcularW(Bcr, Bsr, v, Dx, R);

		calcularBxt(nb, v, Xm, Dx, R, W, t, listBxt);
		List<BigDecimal> valoresCwxt = calcularCwXT(nb, t, listBxt);
		calcularConcentracaoAbaixoFonte(nb, Owv, Oav, pb, koc, foc, kh, valoresCwxt);
		calcularConcentracaoChegaNaAguaSubterranea(nb, q, valoresCwxt, xSource,	ySource);
		calcularPerdas(nb, si, cti, mwTPH, mwi, cTPH, selfi, q, Owv, Oav, pb, koc, foc, kh, pbGL, cw, Defft, Lw, Ld, Bsr, xSource, ySource);
	}



	private static BigDecimal calcularParametroDeDistribuicaoDeTamanhoDeProZonaVadosa(BigDecimal Nv) {
		double powYv = Math.pow(0.5, Nv.divide(Nv.subtract(new BigDecimal("1")), MathContext.DECIMAL128).doubleValue());   //(Nv-1));
		BigDecimal Yv = new BigDecimal("3").add(new BigDecimal("2").divide((Nv.subtract(new BigDecimal("1"))).multiply(new BigDecimal(1).subtract(new BigDecimal(powYv))), MathContext.DECIMAL128));
		return Yv;
	}

	private static BigDecimal calcularConcentracaoInicialFaseDisolvida(BigDecimal cti, BigDecimal Owv, BigDecimal Oav, BigDecimal koc, BigDecimal foc,BigDecimal kh, BigDecimal pbGL) {
		BigDecimal cw = cti.multiply(pbGL).divide(pbGL.multiply(koc).multiply(foc).add(Owv).add(Oav.multiply(kh)), MathContext.DECIMAL128).multiply(new BigDecimal(Math.pow(10, 3))).setScale(30, RoundingMode.HALF_EVEN);
		return cw;
	}

	private static BigDecimal calcularConcentracaoInicialDaFonte(BigDecimal si, BigDecimal selfi, BigDecimal cw) {
		List<BigDecimal> listValues = new ArrayList<BigDecimal>();
		listValues.add(0, cw);
		listValues.add(1, si);
		listValues.add(2, selfi);
		BigDecimal maior = listValues.get(0);
		for (int i = 0; i < listValues.size(); i++) {
			if (maior.compareTo(listValues.get(i)) >= 0) {
				maior = listValues.get(i);
			} 
		}
		BigDecimal cwO = maior;
		return cwO;
	}

	private static void calcularConcentracaoDeVaporInicial(BigDecimal kh,BigDecimal cwO) {
		BigDecimal cvO = kh.multiply(cwO).setScale(30, RoundingMode.HALF_EVEN);
	}

	private static BigDecimal calcularPermeabilidadeRelativa(BigDecimal q) {
		BigDecimal ksu;
		ksu = new BigDecimal("5").multiply(new BigDecimal("100")).setScale(30, RoundingMode.HALF_EVEN);
		BigDecimal Krl = BigDecimal.ZERO;
		Krl = q.divide(ksu).setScale(30, RoundingMode.HALF_EVEN);
		return Krl;
	}

	private static BigDecimal calcularParametroDeDistribuicaoTamanhoDoPoroZonaLens(BigDecimal Nl) {
		double powY = (Math.pow(0.5, Nl.divide(Nl.subtract(new BigDecimal(1)), MathContext.DECIMAL128).doubleValue()));
		BigDecimal yL = new BigDecimal(3).add(new BigDecimal(2).divide((Nl.subtract(new BigDecimal(1)).multiply(new BigDecimal(1).subtract(new BigDecimal(powY)))), MathContext.DECIMAL128));
		return yL;
	}

	private static BigDecimal calcularPorosidadeDaAgua(BigDecimal Krl, BigDecimal Orl, BigDecimal Otl, BigDecimal yL) {
		double powOwl = Math.pow(Krl.doubleValue(), new BigDecimal("1").divide(yL, MathContext.DECIMAL128).doubleValue());
		BigDecimal Owl = new BigDecimal(powOwl).multiply(Otl.subtract(Orl)).add(Orl).setScale(30, RoundingMode.HALF_EVEN);
		return Owl;
	}

	private static BigDecimal calcularPorosidadeDoArZonaLens(BigDecimal Otl, BigDecimal Owl) {
		BigDecimal Oal = Otl.subtract(Owl).setScale(30, RoundingMode.HALF_EVEN);
		return Oal;
	}

	private static BigDecimal calcularCoeficienteDeDifusaoEfetivaZonaVadosa(BigDecimal Otv, BigDecimal Owv, BigDecimal Oav, BigDecimal kh, BigDecimal Dw, BigDecimal Da) {
		double powOav = Math.pow(Oav.doubleValue(), 3.333333333);
		double powOtv = Math.pow(Otv.doubleValue(), 2d);
		double powOwv = Math.pow(Owv.doubleValue(), 3.333333333);
		
		BigDecimal Deffv = (Da.multiply(new BigDecimal(powOav).divide(new BigDecimal(powOtv), MathContext.DECIMAL128)).add(Dw.multiply(new BigDecimal("1").divide(kh, MathContext.DECIMAL128)).multiply(new BigDecimal(powOwv).divide(new BigDecimal(powOtv), MathContext.DECIMAL128))).setScale(30, RoundingMode.HALF_EVEN));
		return Deffv;
	}

	private static BigDecimal calcularCoeficienteDeDifusaoEfetivaZonaLens(BigDecimal kh, BigDecimal Otl, BigDecimal Owl, BigDecimal Oal, BigDecimal Dw, BigDecimal Da) {
		double powOwl;
		double powOal = Math.pow(Oal.doubleValue(), 3.333333333);
		double powOtl = Math.pow(Otl.doubleValue(), 2d);
		powOwl = Math.pow(Owl.doubleValue(), 3.333333333);
		
		BigDecimal Deffl = (Da.multiply(new BigDecimal(powOal).divide(new BigDecimal(powOtl), MathContext.DECIMAL128)).add(Dw.multiply(new BigDecimal("1").divide(kh, MathContext.DECIMAL128)).multiply(new BigDecimal(powOwl).divide(new BigDecimal(powOtl), MathContext.DECIMAL128))).setScale(30, RoundingMode.HALF_EVEN));
		return Deffl;
	}

	private static BigDecimal calcularCoeficienteDeDifusaoEfetivaTotal(BigDecimal Deffv, BigDecimal Deffl, BigDecimal L1, BigDecimal Lv,BigDecimal Lv1) {
		BigDecimal Defft = (Lv.add(L1)).divide(L1.divide(Deffl).add(Lv1.divide(Deffv, MathContext.DECIMAL128)), MathContext.DECIMAL128);
		return Defft;
	}
	
	private static BigDecimal conveterSiParaCM3(BigDecimal si) {
		BigDecimal siGcm3 = si.multiply(new BigDecimal(Math.pow(10d, -6d))).setScale(9, RoundingMode.UP);
		return siGcm3;
	}

	private static BigDecimal calcularDistanciaSuperficieDoSoloCentroDaFonte(BigDecimal Lw, BigDecimal Psource) {
		BigDecimal Ld;
		Ld = Psource.add(Lw.divide(new BigDecimal("2"))).multiply(new BigDecimal("100"));
		return Ld;
	}

	private static BigDecimal converterDefftParaCm2PorDia(BigDecimal Defft) {
		BigDecimal defftCm2D = Defft.multiply(new BigDecimal("24")).multiply(new BigDecimal("60")).multiply(new BigDecimal("60"));
		return defftCm2D;
	}

	private static BigDecimal converterLwParaCm(BigDecimal Lw) {
		BigDecimal LwCM = Lw.multiply(new BigDecimal("100"));
		return LwCM;
	}

	private static BigDecimal calcularQuedaDaFontePorPercolacaoComFaseResidual(BigDecimal mwTPH, BigDecimal mwi, BigDecimal cTPH, BigDecimal q, BigDecimal pb, BigDecimal siGcm3, BigDecimal LwCM) {
		BigDecimal Bcrp;
		Bcrp = (q.multiply(mwTPH).multiply(siGcm3)).divide(pb.multiply(LwCM).multiply(cTPH).multiply(mwi), MathContext.DECIMAL128);
		return Bcrp;
	}

	private static BigDecimal calcularQuedaDaFonteSemFaseResidual(BigDecimal Bsrp, BigDecimal Bsrv) {
		BigDecimal Bsr;
		Bsr = Bsrp.add(Bsrv);
		return Bsr;
	}

	private static BigDecimal calcularQuedaDaFontePorVolatizacaoComFaseResidual(BigDecimal mwTPH, BigDecimal mwi, BigDecimal cTPH, BigDecimal pb,BigDecimal kh, BigDecimal siGcm3, BigDecimal Ld, BigDecimal defftCm2D, BigDecimal LwCM) {
		BigDecimal Bcrv;
		Bcrv = (defftCm2D.multiply(kh).multiply(mwTPH).multiply(siGcm3)).divide(pb.multiply(LwCM).multiply(Ld).multiply(mwi).multiply(cTPH), MathContext.DECIMAL128);
		System.out.println("Valor de Bcrv: " + Bcrv);
		return Bcrv;
	}

	private static BigDecimal calcularQuedaDaFontePorPercolacaoSemFaseResidual(NumberFormat nb, BigDecimal q, BigDecimal Owv, BigDecimal Oav, BigDecimal KOC, BigDecimal pb, BigDecimal foc, BigDecimal kh,	BigDecimal LwCM) {
		BigDecimal Bsrp;
		Bsrp  = q.divide((pb.multiply(KOC).multiply(foc)).add(Owv.add(Oav.multiply(kh))).multiply(LwCM), MathContext.DECIMAL128);
		return Bsrp;
	}

	private static BigDecimal calcularQuedaDaFontePorVolatizacaoSemFaseResidual(NumberFormat nb, BigDecimal Owv, BigDecimal Oav, BigDecimal KOC, BigDecimal pb, BigDecimal foc, BigDecimal kh, BigDecimal Ld, BigDecimal defftCm2D, BigDecimal LwCM) {
		BigDecimal Bsrv;
		Bsrv = (defftCm2D.multiply(kh)).divide(((pb.multiply(KOC).multiply(foc)).add(Owv.add(Oav.multiply(kh))).multiply(LwCM.multiply(Ld))), MathContext.DECIMAL128);
		return Bsrv;
	}

	private static BigDecimal calcularW(BigDecimal Bcr, BigDecimal Bsr, BigDecimal v, BigDecimal Dx, BigDecimal R) {
		BigDecimal W = BigDecimal.ZERO; 
		BigDecimal U = new BigDecimal("0.0034");
	
		if (Bsr.compareTo(new BigDecimal("0")) >= 0) {
			W = v.multiply(new BigDecimal(Math.pow((new BigDecimal(1)).add((new BigDecimal(4).multiply(Dx)).divide(new BigDecimal(Math.pow(v.doubleValue(), 2)), MathContext.DECIMAL128).multiply(U.subtract(R.multiply(Bsr)))).doubleValue(), 1d / 2d)));
		}

		if (Bcr.compareTo(new BigDecimal("0")) >= 0) { //TODO: Sempre entrando aqui..
//			W = v.multiply(new BigDecimal(Math.pow(
//					(new BigDecimal(1)).add(
//							(new BigDecimal(4).multiply(Dx))
//									.divide(new BigDecimal(Math.pow(
//											v.doubleValue(), 2)),
//											MathContext.DECIMAL128).multiply(
//											U.subtract(R.multiply(Bcr))))
//							.doubleValue(), 1d / 2d)));
//			System.out.println("Valor de W com Bcr () : " + W.setScale(9, RoundingMode.HALF_EVEN));
//			W = v * Math.sqrt(1 + ((4 * Dx)/Math.pow(v, 2)) * (U - R.doubleValue() *  Bcr ));
		}
		return W;
	}

	private static BigDecimal calcularFatorDeRetardacao(BigDecimal Owv,	BigDecimal KOC, BigDecimal pb, BigDecimal foc) {
		BigDecimal R = new BigDecimal("1").add((pb.multiply(foc).multiply(KOC)).divide(Owv, MathContext.DECIMAL128)).setScale(20, RoundingMode.HALF_EVEN);
		return R;
	}

	private static BigDecimal calcularDoCoeficienteDeDispersao(BigDecimal Xm) {
		BigDecimal aLx = BigDecimal.ZERO;
		if (Xm.compareTo(new BigDecimal("2")) >= 0) {
			aLx = new BigDecimal(Math.exp(0.584 * Math.log(Xm.doubleValue()) - 2.727));
		} else  {
			aLx = new BigDecimal(Math.exp((-1)*4.933+3.811*Math.log(Xm.doubleValue())));
		}
		return aLx;
	}

	private static BigDecimal calcularDistanciaDaFonteAteLencolFreatico(BigDecimal Lv, BigDecimal Lw, BigDecimal Psource) {
		BigDecimal Xm = Lv.subtract(Lw).subtract(Psource).setScale(30, RoundingMode.HALF_EVEN);
		return Xm;
	}

	private static BigDecimal calcularVelocidadeDePercolacao(BigDecimal q, BigDecimal Owv) {
		BigDecimal v = q.divide(Owv, MathContext.DECIMAL128);
		return v;
	}

	private static void calcularBxt(NumberFormat nb, BigDecimal v, BigDecimal Xm, BigDecimal Dx, BigDecimal R, BigDecimal W, BigDecimal t, List<BigDecimal> listBxt) {
		BigDecimal x = BigDecimal.ZERO;

		/**
		 * Convertendo valor de v para m/d
		 */
		BigDecimal Vmd = v.multiply(new BigDecimal(Math.pow(10d, -2d))).setScale(30, RoundingMode.HALF_EVEN); 
		
		/**
		 * convertendo valor de W para m/d
		 */
		BigDecimal Wmd = W.divide(new BigDecimal("100"), MathContext.DECIMAL128).setScale(30, RoundingMode.HALF_EVEN);

		/**
		 * Convertendo DX para m³/d
		*/
		BigDecimal Dxm3d = Dx.divide(new BigDecimal("10000")).setScale(30, RoundingMode.HALF_EVEN);
		Integer valorTermo = 1;
		
		while (x.compareTo(Xm) <= 0) {
			BigDecimal CapilaryZone = VADOSE_ZONE.subtract(DEPTH_TO_TOP).subtract(THICKNESS_OF_SOURCE);
			BigDecimal valueX = x.add(new BigDecimal("0.8"));
			x = valorTermo == 1 ? BigDecimal.ZERO : valueX.compareTo(Xm)>0 ? CapilaryZone : valueX;
			
			/**
	 		* Cálculo do Primeiro Termo
			*/
			BigDecimal primeiroTermo = (Vmd.subtract(Wmd).multiply(x)).divide(Dxm3d.multiply(new BigDecimal(2d)), 20, RoundingMode.HALF_UP);
				
			/**
			 * Cálculo do segundoTermo
			 */
			BigDecimal dividendo = R.multiply(x).subtract(Wmd.multiply(t));
			BigDecimal dois = new BigDecimal("2").setScale(20);
			BigDecimal divisor = new BigDecimal(Math.sqrt(Dxm3d.multiply(R).multiply(t).doubleValue())).multiply(dois);
			BigDecimal segundoTermo = dividendo.divide(divisor, MathContext.DECIMAL128).setScale(30, RoundingMode.HALF_EVEN);
				
			/**
			 * Cálculo do terceiroTermo
			 */
			BigDecimal terceiroTermo =  (Vmd.add(Wmd).multiply(x)).divide(Dxm3d.multiply(new BigDecimal(2d)),20,RoundingMode.HALF_EVEN);
	
			/**
			 * Cálculo do quartoTermo
			 */
			BigDecimal quartoTermo = (R.multiply(x).add(Wmd.multiply(t)).divide((bigSqrt(Dxm3d.multiply(R).multiply(t)).multiply(dois)), MathContext.DECIMAL128).setScale(30, RoundingMode.HALF_EVEN));
			BigDecimal divisao = new BigDecimal(1d/2d);
			BigDecimal BxtFinal = (divisao.multiply(new BigDecimal(Math.exp(primeiroTermo.doubleValue()))).multiply(new BigDecimal(Erf.erfc(segundoTermo.doubleValue()))).add((divisao).multiply(new BigDecimal(Math.exp(terceiroTermo.doubleValue()))).multiply(new BigDecimal(Erf.erfc(quartoTermo.doubleValue())))));
			listBxt.add(BxtFinal);
			if(valueX.compareTo(Xm) > 0){
				break;
			}
			valorTermo++;
		}
	}

	private static List<BigDecimal> calcularCwXT(NumberFormat nb, BigDecimal t,	List<BigDecimal> listBxt) {
		BigDecimal CWO = new BigDecimal("34.104417354216200000000000000000");
		BigDecimal Beta = new BigDecimal("0.000791788970058695000000000000");
		List<BigDecimal> valoresCwxt = new ArrayList<BigDecimal>();
		
		for (BigDecimal btx : listBxt) {
			BigDecimal Cwxt = CWO.multiply(new BigDecimal(Math.exp(-Beta.multiply(t).doubleValue())).multiply(btx));
			valoresCwxt.add(Cwxt);
		}
		return valoresCwxt;
	}

	private static void calcularConcentracaoAbaixoFonte(NumberFormat nb, BigDecimal Owv, BigDecimal Oav, BigDecimal pb, BigDecimal koc, BigDecimal foc, BigDecimal kh, List<BigDecimal> valoresCwxt) {
		ArrayList<BigDecimal> valoresCt = new ArrayList<BigDecimal>();
		Integer iCt = 1;
		
		for (BigDecimal cwxt : valoresCwxt) {
			BigDecimal Ct = (cwxt.multiply(((pb.multiply(new BigDecimal(Math.pow(10d, 3d)))).multiply(koc).multiply(foc)).add(Owv).add(Oav.multiply(kh))).divide(pb.multiply(new BigDecimal(Math.pow(10d, 3d))), MathContext.DECIMAL128)).multiply(new BigDecimal(Math.pow(10d, 3d)));
			valoresCt.add(Ct);
			iCt++;
		}
	}

	private static void calcularConcentracaoChegaNaAguaSubterranea(NumberFormat nb, BigDecimal q, List<BigDecimal> valoresCwxt,BigDecimal xSource, BigDecimal ySource) {
		BigDecimal Mass = BigDecimal.ZERO;
		BigDecimal C = valoresCwxt.get(valoresCwxt.size() - 1); // é o ultimo valor de Cwxt
		Mass = q.multiply(C).divide(xSource.multiply(ySource), MathContext.DECIMAL128).setScale(30, RoundingMode.HALF_EVEN);
	}

	private static void calcularPerdas(NumberFormat nb, BigDecimal si, BigDecimal cti, BigDecimal mwTPH, BigDecimal mwi, BigDecimal cTPH, BigDecimal selfi, BigDecimal q, BigDecimal Owv, BigDecimal Oav, BigDecimal pb, BigDecimal koc, BigDecimal foc, BigDecimal kh, BigDecimal pbGL, BigDecimal cw, BigDecimal Defft, BigDecimal Lw, BigDecimal Ld, BigDecimal Bsr, BigDecimal xSource,	BigDecimal ySource) {
		
		BigDecimal t;
		System.out.println("");
		System.out.println("================================================");
		System.out.println("================ PERDAS ========================");
		System.out.println("================================================");
		System.out.println("");
		
		List<BigDecimal> totaisDeMassaNafonte = new ArrayList<BigDecimal>();
		List<BigDecimal> totaisDeCTGG = new ArrayList<BigDecimal>();
		List<BigDecimal> totaisDeCwGcm3 = new ArrayList<BigDecimal>();
		List<BigDecimal> totaisDeCvGcm3 = new ArrayList<BigDecimal>();
		BigDecimal CtGG = BigDecimal.ZERO;
		BigDecimal totalDeMassaNaFonteIndice0 = BigDecimal.ZERO;  //KG
		t = new BigDecimal("0");

		BigDecimal DefftCm2d = converterDefftParaCm2PorDia(Defft);
		
		for (int ty = 0; ty <= 10; ty++) {

			if (ty == 0 && t.equals(new BigDecimal("0"))) {
				totalDeMassaNaFonteIndice0 = cti.multiply(pb).multiply(xSource).multiply(ySource).multiply(Lw.multiply(new BigDecimal("100"))).multiply(new BigDecimal(Math.pow(10d, -3d))).setScale(30, RoundingMode.HALF_EVEN);
				CtGG = totalDeMassaNaFonteIndice0.divide(pb.multiply(xSource).multiply(ySource).multiply(Lw.multiply(new BigDecimal("100"))), MathContext.DECIMAL128).multiply(new BigDecimal(Math.pow(10d, 3d))).setScale(30, RoundingMode.HALF_EVEN);
			}
			
			BigDecimal tmf = totalDeMassaNaFonteIndice0.multiply(new BigDecimal(Math.exp(-Bsr.multiply(t).doubleValue())));
			System.out.println("Valor de TMF (" + ty + "): ====> " + nb.format(tmf));
			totaisDeMassaNafonte.add(tmf);
			
			//********************************** Ct(g/g) ***************************************//
			CtGG = totaisDeMassaNafonte.get(ty).divide(pb.multiply(xSource).multiply(ySource).multiply(Lw.multiply(new BigDecimal("100"))), MathContext.DECIMAL128).multiply(new BigDecimal(Math.pow(10d, 3d))).setScale(30, RoundingMode.HALF_EVEN);
			System.out.println("Valor de Ct(g/g): " + nb.format(CtGG));
			totaisDeCTGG.add(CtGG);
			
			//************************************ CW(g/cm³) ***********************************//
			BigDecimal CwGcm3 = BigDecimal.ZERO;
			if (cw.compareTo(selfi) >= 0 ) {
				CwGcm3 = totaisDeCTGG.get(ty).multiply(mwTPH).divide(mwi.multiply(cTPH).multiply(si), MathContext.DECIMAL128).multiply(new BigDecimal(Math.pow(10d, -6d))); //TODO: Verificar com a natalia 
				System.out.println("Valor de CwGcm3: " + nb.format(CwGcm3));
				totaisDeCwGcm3.add(CwGcm3);
			} else {
				CwGcm3 = CtGG.multiply(pbGL).divide((pbGL.multiply(koc).multiply(foc)).add(Owv.add(Oav.multiply(kh))), MathContext.DECIMAL128).multiply(new BigDecimal(Math.pow(10d, -3d)));
				System.out.println("Valor de CwGcm3: " + nb.format(CwGcm3));
				totaisDeCwGcm3.add(CwGcm3);
			}

			//************************************ Cv(g/cm³) ************************************//
			BigDecimal CvGcm3 = BigDecimal.ZERO;
			CvGcm3 = kh.multiply(totaisDeCwGcm3.get(ty));
			totaisDeCvGcm3.add(CvGcm3);
			System.out.println("Valor de CvGcm3: " + nb.format(CvGcm3));
			
			/**
			 * Liq. Mass Lost From Source (kg)
			 */
			BigDecimal ASource = xSource.multiply(ySource);
			BigDecimal MassLot = BigDecimal.ZERO;
			MassLot = q.multiply(totaisDeCwGcm3.get(ty)).multiply(ASource.multiply(t)).multiply(new BigDecimal(Math.pow(10d, -3d)));
			System.out.println("Valor de MassLot: " + nb.format(MassLot));
			
			/**
			 * Volatization loss (kg)
			 */
			BigDecimal volationLoss = BigDecimal.ZERO;
			volationLoss = ((DefftCm2d.multiply(totaisDeCvGcm3.get(ty).divide(Ld))).multiply(t).multiply(ASource)).multiply(new BigDecimal(Math.pow(10d, -3d)));
			System.out.println("Valor de voltationLoss: " + nb.format(volationLoss));
			
			System.out.println("\n=========================================");
			t = t.add(new BigDecimal("365"));
		}
	}
	
	/**
	 * Private utility method used to compute the square root of a BigDecimal.
	 * 
	 * @author Luciano Culacciatti 
	 * @url http://www.codeproject.com/Tips/257031/Implementing-SqrtRoot-in-BigDecimal
	 */
	private static BigDecimal sqrtNewtonRaphson  (BigDecimal c, BigDecimal xn, BigDecimal precision){
	    BigDecimal fx = xn.pow(2).add(c.negate());
	    BigDecimal fpx = xn.multiply(new BigDecimal(2));
	    BigDecimal xn1 = fx.divide(fpx,2*SQRT_DIG.intValue(),RoundingMode.HALF_DOWN);
	    xn1 = xn.add(xn1.negate());
	    BigDecimal currentSquare = xn1.pow(2);
	    BigDecimal currentPrecision = currentSquare.subtract(c);
	    currentPrecision = currentPrecision.abs();
	    if (currentPrecision.compareTo(precision) <= -1){
	        return xn1;
	    }
	    return sqrtNewtonRaphson(c, xn1, precision);
	}

	/**
	 * Uses Newton Raphson to compute the square root of a BigDecimal.
	 * 
	 * @author Luciano Culacciatti 
	 * @url http://www.codeproject.com/Tips/257031/Implementing-SqrtRoot-in-BigDecimal
	 */
	public static BigDecimal bigSqrt(BigDecimal c){
	    return sqrtNewtonRaphson(c,new BigDecimal(1),new BigDecimal(1).divide(SQRT_PRE));
	}
}