Vadose Zone 1.1

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.HashMap;
import java.util.List;
import java.util.Map;

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;
import br.com.sed.seris.domain.simulation.SiteCharacterizationSoil;
import br.com.sed.seris.domain.simulation.VadoseZoneResultWrapper;

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");
	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 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 S = 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 MW = 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(MW, MathContext.DECIMAL128)).multiply(S).setScale(30, RoundingMode.HALF_EVEN);
		
		System.out.println("Valor de Selfi: " + nb.format(selfi));
		
		// 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 KsuVsas = new BigDecimal(soil.getHydraulicConductivity()).setScale(30, RoundingMode.HALF_EVEN);
		BigDecimal KsuV = 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(KsuV) >= 0 ) {
			KsuV = q;
		}
		
		krv = q.divide(KsuV).setScale(30, RoundingMode.HALF_EVEN);
		System.out.println("Valor de Selfi: " + nb.format(krv));
		
		BigDecimal Yv = calcularParametroDeDistribuicaoDeTamanhoDeProZonaVadosa(Nv);
		System.out.println("Valor de Selfi: " + nb.format(Yv));
		
		BigDecimal cm3 = new BigDecimal("1");
		BigDecimal OwV = calcularPorosidadeDaAgua(krv, OrV, OtV, Yv);
		System.out.println("Valor de Owv: " + nb.format(OwV));
		
		// 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);
		System.out.println("Valor de Kh: " + nb.format(Kh));
		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);
		System.out.println("Valor de Owv: " + nb.format(cw));
		
		//PASSO 1.2.3 : Concentração inicial na fonte (Cw0)
		BigDecimal cwO = calcularConcentracaoInicialDaFonte(S, selfi, cw);
		System.out.println("Valor de cw0" + nb.format(cwO));
		
		// Passo 2 : Cálculo da concentração de vapor Inicial (Cv0)
		calcularConcentracaoDeVaporInicial(Kh, cwO);
				
		// PASSO 3 :  CÁLCULO DE BETA (β)
		BigDecimal Krl = calcularPermeabilidadeRelativa(q);
		System.out.println("Valor de KRL:" + nb.format(Krl));
				
		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);
		System.out.println("Valor de yL" + nb.format(yL));
		
		BigDecimal OwL = calcularPorosidadeDaAgua(Krl, OrL, OtL, yL);
		System.out.println("OWL" + nb.format(OwL));
		
		BigDecimal OaL = calcularPorosidadeDoArZonaLens(OtL, OwL);
		System.out.println("OAL" + nb.format(OaL));
		
		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);
		System.out.println("DEFFV" + nb.format(Deffv));
		
		BigDecimal Deffl = calcularCoeficienteDeDifusaoEfetivaZonaLens(Kh, OtL,	OwL, OaL, Dw, Da);
		System.out.println("DEFFL" + nb.format(Deffl));

		BigDecimal Psource = new BigDecimal("2"); //TODO:  entrada do usuário --> depth to top
		BigDecimal hL = new BigDecimal("0"); //m //TODO: entrada do usuário --> Thickness of lens zone
		BigDecimal hV = new BigDecimal("10"); //TODO: entrada do usuário --> Thickness of vadose zone
		BigDecimal hS = Psource.subtract(hL);
		
		//BigDecimal Lv1 = hV.subtract(hL);
		BigDecimal Defft = calcularCoeficienteDeDifusaoEfetivaTotal(Deffv, Deffl, hL, hS);
		
		System.out.println("Valor de Defft: " + nb.format(Defft));
	    
		BigDecimal Lw = new BigDecimal("1"); //TODO: declaracao
		BigDecimal S_gcm3 = conveterSiParaCM3(S);

		BigDecimal Ld = Psource.add(Lw.divide(new BigDecimal("2")));
		BigDecimal Defft_cm2_d = converterDefftParaCm2PorDia(Defft);
		BigDecimal Lw_cm = converterLwParaCm(Lw);
		Ld = calcularDistanciaSuperficieDoSoloCentroDaFonte(Lw, 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, MW,	Ctph, q, pb, S_gcm3, Lw_cm);
			Bcrv = calcularQuedaDaFontePorVolatizacaoComFaseResidual(MWtph, MW, Ctph, pb, Kh, S_gcm3, Ld, Defft_cm2_d, Lw_cm);
			Bcr = calcularQuedaDaFonteComFaseResidual(Bcrp, Bcrv);
		} else {
			Bsrp = calcularQuedaDaFontePorPercolacaoSemFaseResidual(nb, q, OwV, OaV, Koc, pb, Foc, Kh, Lw_cm);
			Bsrv = calcularQuedaDaFontePorVolatizacaoSemFaseResidual(nb, OwV, OaV, Koc, pb, Foc, Kh, Ld, Defft_cm2_d, Lw_cm);
			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(hV, Lw, Psource);
		BigDecimal Xm1 = new BigDecimal("1");
		BigDecimal XmO = Xm.divide(Xm1);
		BigDecimal aL = calcularDoCoeficienteDeDispersao(Xm);
		BigDecimal Dx = aL.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);
		
		System.out.println("\n============================================");
		System.out.println("============ Valores de CWXT =================");
		System.out.println("==============================================");
		
		System.out.println("Tamanho" + valoresCwxt.size());
		
		for (int i = 0; i < valoresCwxt.size(); i++) {
			System.out.println(nb.format(valoresCwxt.get(i)));
		}
		
		calcularConcentracaoAbaixoFonte(nb, OwV, OaV, pb, koc, Foc, Kh, valoresCwxt);
		List<BigDecimal> valoresMass = calcularConcentracaoChegaNaAguaSubterranea(nb, q, valoresCwxt, xSource,	ySource);
		
		System.out.println("\n============================================");
		System.out.println("============ Valores de Mass =================");
		System.out.println("==============================================");
		
		for (int i = 0; i < valoresMass.size(); i++) {
			System.out.println(nb.format(valoresMass.get(i)));
		}
		
		calcularPerdas(nb, S, Cti, MWtph, MW, 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 hL, BigDecimal hV) {
		BigDecimal Defft = (hV.add(hL)).divide(hL.divide(Deffl, MathContext.DECIMAL128).add(hV.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 calcularQuedaDaFonteComFaseResidual(BigDecimal Bcrp, BigDecimal Bcrv) {
		BigDecimal Bsr;
		Bsr = Bcrp.add(Bcrv);
		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 hV, BigDecimal Lw, BigDecimal Psource) {
		BigDecimal Xm = hV.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) {
		
		
		System.out.println("\n===============================================");
		System.out.println("=============== B (x,t ) ======================");
		System.out.println("===============================================\n");
		
		
		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;
		
		for (int i = 0; i < 10; i ++) {
			System.out.println("================= " + i);
			
			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++;
				System.out.println("Valor de BxtFinal: " + nb.format(BxtFinal));
			}
	
			x = new BigDecimal("0.8");
			t = t.add(new BigDecimal("365"));
		
		}
		
	}

	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 List<BigDecimal> calcularConcentracaoChegaNaAguaSubterranea(NumberFormat nb, BigDecimal q, List<BigDecimal> valoresCwxt,BigDecimal xSource, BigDecimal ySource) {
		List<BigDecimal> valoresMass = new ArrayList<BigDecimal>();
		
		for (int i = 0; i < valoresCwxt.size(); i++) {
			BigDecimal Mass = BigDecimal.ZERO;
			Mass = q.multiply(valoresCwxt.get(i)).divide(xSource.multiply(ySource), MathContext.DECIMAL128).setScale(30, RoundingMode.HALF_EVEN);
			valoresMass.add(Mass);
		}
		
		return valoresMass;
	}

	private static List<VadoseZoneResultWrapper> 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 yearSimulation = BigDecimal.ZERO;
		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>();
		List<VadoseZoneResultWrapper> vadoseZoneWrapper = new ArrayList<VadoseZoneResultWrapper>();

		BigDecimal CtGG = BigDecimal.ZERO;
		BigDecimal totalDeMassaNaFonteIndice0 = BigDecimal.ZERO; 
		BigDecimal DefftCm2d = converterDefftParaCm2PorDia(Defft);
		
		for (int ty = 0; ty <= 10; ty++) { //TODO: Time your simulation

			VadoseZoneResultWrapper vadoseZoneResultWrapper = new VadoseZoneResultWrapper();
			vadoseZoneResultWrapper.setTimeSimulation((long)ty);
			
			if (ty == 0 && yearSimulation.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(yearSimulation).doubleValue())));
			System.out.println("Valor de TMF (" + ty + "): ====> " + nb.format(tmf));
			totaisDeMassaNafonte.add(tmf);
			vadoseZoneResultWrapper.setTmf(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);
			vadoseZoneResultWrapper.setCtg(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);
			}
			vadoseZoneResultWrapper.setCwxt(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));
			vadoseZoneResultWrapper.setCv(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(yearSimulation)).multiply(new BigDecimal(Math.pow(10d, -3d)));
			System.out.println("Valor de MassLot: " + nb.format(MassLot));
			vadoseZoneResultWrapper.setMass(MassLot);
			
			/**
			 * Volatization loss (kg)
			 */
			BigDecimal volationLoss = BigDecimal.ZERO;
			volationLoss = ((DefftCm2d.multiply(totaisDeCvGcm3.get(ty).divide(Ld))).multiply(yearSimulation).multiply(ASource)).multiply(new BigDecimal(Math.pow(10d, -3d)));
			System.out.println("Valor de voltationLoss: " + nb.format(volationLoss));
			vadoseZoneResultWrapper.setVolationLoss(volationLoss);
			
			vadoseZoneWrapper.add(vadoseZoneResultWrapper);
			vadoseZoneResultWrapper = new VadoseZoneResultWrapper();
			System.out.println("\n=========================================");
			yearSimulation = yearSimulation.add(new BigDecimal("365"));
		}
		
		return vadoseZoneWrapper;
	}
	
	/**
	 * 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));
	}
}