/* MA-CME Notes:
Equivalent input line:
rules file=NTS-LV-Rules.java Qf=1 minthick=0.01 minVs=0.35 minVp=0.606 maxv=7.8 vpovs=1.732051 zgeotech=0.03 rockgeotech=0.76 soilgeotech=0.5
Southern Nevada regional model, with earthquake-location Vp(z) for rock and 
Jachens, Moring, Blakely den(z) for both volcanic and sedimentary basins.
Q relations from Olsen and Day for LA models.
*/

import java.lang.Math;
import java.util.*;

public class Rules
	{
	/* Version 5 additions */
	static String rulesline="";
	static String rulessource="";
	static String rulesdir="";
	/* Q value considered by E3D to mean purely linear wave propagation */
	static float noatten=2000F;
	static float Qf=1F; /* Central frequency of attenuation curve */
	/* Minimum sediment thickness to be considered in a basin, km */
	static float minthick=0.01F;
	/* Minimum shear velocity for any grid node, km/s */
	static float minVs=0.35F;
	/* Minimum P velocity for any grid node, km/s */
	static float minVp=0.606F;
	/* Maximum velocity for any grid node, km/s */
	static float maxv=7.8F;
	/* Version 3 values below */
	/* Default Vp/Vs ratio */
	static float vpovs=1.732051F; /* Poisson solid */
	static float zgeotech=0.03F; /* km, default geotech velocities to 30 m depth */
	static float rockgeotech=0.76F;
		/* km/s, default 30-m shear velocity in rock, NEHRP B-C Boundary */
	static float soilgeotech=0.5F;
		/* default Vs30 in most of Las Vegas Valley */
	/* Classic relation from Gardner, G.H.F., L.W. Gardner, and A.R. Gregory, 1974, Formation
		velocity and density; the diagnostic basics for stratigraphic traps, Geophysics, 39,
		770Ð780. */
	static float vpfromden(float den) /* km/s from g/cc */
		{
		double qv;
		qv = den/0.23;
		return (float)(qv*qv*qv*qv*0.3048/1000);
		}
	static float vsfromvp(float vp)
		{
		return vp/vpovs;
		}
	static float vpfromvs(float vs)
		{
		return vs*vpovs;
		}
	static float denfromvp(float vp) /* g/cc from km/s */
		{
		if (vp < minVp)
			vp = minVp;
		return (float)(0.23*Math.pow(vp*1000/0.3048, 0.25));
		}

	static float bedvp(float z) /* WGB earthquake location model */
		{
		float v0=5.5F; /* km/s, revised for WGB from 5.0 */
		float v1=6F;
		float z1=5F; /* km */
		float v2=7.8F;
		float z2=35F;
		if (z <= z1)
			return v0 + z*(v1 - v0)/z1;
		else if (z <= z2)
			return v1;
		else
			return v2;
		}

	/* New Version 5 Rules */
	/* Implementing attenuation for Version 5: Olsen et al. (2003);
	   Day (1998) coarse-grained viscoelasticity fitting LA Basin data */
	static float bedQs(float vs)
		{
		return Rules.noatten;
		}
	static float bedQp(float Qs)
		{
		return Rules.noatten;
		}
	static float basinQs(float vs)
		{
		if (vs < 1.5F)
			return 20F*vs;
		return 100F*vs;
		}
	static float basinQp(float Qs)
		{
		return 1.5F*Qs;
		}

	/* Bedrock velocity finder disregards lat-lon and only uses Rules on z */
	static ElastProps getBedrockProps(Grid grid, double lat, double lon,
		float z, float rockvel)
		{
		float den, vp, vs, Qp, Qs;
		Qs = Qp = -1F;
		if (z < grid.dh || z < zgeotech)
			return Rules.getRockGeotechProps(Rules.zgeotech,
				rockvel, grid.dh, grid.atten, z);
		vp = Rules.bedvp(z);
		vs = Rules.vsfromvp(vp);
		den = Rules.denfromvp(vp);
		/* Implement Q from Benites & Olsen, 2005 (WnLH & LA) */
		if (grid.atten)
			{
			Qs = Rules.bedQs(vs);
			Qp = Rules.bedQp(Qs);
			}
		/* Implement minimum velocities */
		if (vp < Rules.minVp)
			vp = Rules.minVp;
		if (vs < Rules.minVs)
			vs = Rules.minVs;
		if (vs > Rules.maxv)
			vs = Rules.maxv;
		if (vp > Rules.maxv)
			vp = Rules.maxv;
		if (grid.atten)
			return new ElastProps(vp, vs, den, Qp, Qs);
		return new ElastProps(vp, vs, den);
		}

	/* Basin velocity finder disregards lat-lon and only uses Rules on z */
	static ElastProps getBasinProps(Grid grid, double lat, double lon,
		float thick, float z, float soilvel, float geol)
		{
		float den, vp, vs, Qp, Qs;
		Qs = Qp = -1F;
		if (z < grid.dh || z < zgeotech)
			return Rules.getSoilGeotechProps(Rules.zgeotech,
			Rules.rockgeotech, soilvel, grid.dh, grid.atten, thick, z);
		if (Rules.isVolcanic(geol))
			{
			den = Rules.volcbasinden(z);
			vp = Rules.vpfromden(den);
			vs = Rules.vsfromvp(vp);
			}
		else
			{
			den = Rules.basinden(z);
			vp = Rules.vpfromden(den);
			vs = Rules.vsfromvp(vp);
			}
		/* Implement Q from Benites & Olsen, 2005 (WnLH & LA) */
		if (grid.atten)
			{
			Qs = Rules.basinQs(vs);
			Qp = Rules.basinQp(Qs);
			}
		/* Implement minimum velocities */
		if (vp < Rules.minVp)
			vp = Rules.minVp;
		if (vs < Rules.minVs)
			vs = Rules.minVs;
		if (vs > Rules.maxv)
			vs = Rules.maxv;
		if (vp > Rules.maxv)
			vp = Rules.maxv;
		if (grid.atten)
			return new ElastProps(vp, vs, den, Qp, Qs);
		return new ElastProps(vp, vs, den);
		}

	static ElastProps getSoilGeotechProps(float dh, boolean atten, float thick, float z)
		{
		return getSoilGeotechProps(zgeotech, rockgeotech, soilgeotech,
			dh, atten, thick, z);
		}
	static ElastProps getSoilGeotechProps(float gtdepth, float rockvel, float soilvel,
		float dh, boolean atten, float thick, float z)
		{
		float den, vp, vs, vg, vb, vr, Qs, Qp;
		Qs = Qp = -1F;
		if (z >= dh && z > gtdepth) /* geotech only affects surface of grid */
			{
			den = Rules.basinden(z);
			vp = Rules.vpfromden(den);
			vs = Rules.vsfromvp(vp);
			}
		else if (dh <= gtdepth) /* use only geotech velocity */
			{
			if (thick < gtdepth)
				vs = rockvel;
			else
				vs = soilvel;
			vp = Rules.vpfromvs(vs);
			den = Rules.denfromvp(vp);
			}
		else /* dh > gtdepth, use slowness-avg surface bedrock w/ geotech */
			{
			if (thick <= gtdepth)
				vg = gtdepth/rockvel;
			else
				vg = gtdepth/soilvel;
			if (thick < dh)
				{
				if (thick <= gtdepth)
					{
					vr = (dh-gtdepth)/(Rules.vsfromvp(
						Rules.bedvp(0F)));
					vb = 0F;
					}
				else
					{
					vr = (dh-thick)/(Rules.vsfromvp(
						Rules.bedvp(0F)));
					vb = (thick-gtdepth)/(Rules.vsfromvp(
					Rules.vpfromden(Rules. basinden(0F))));
					}
				vs = dh/(vg + vb + vr);
				}
			else
				{
				vb = (dh-zgeotech)/(Rules.vsfromvp(
					Rules.vpfromden(Rules. basinden(0F))));
				vs = dh/(vg + vb);
				}
			vp = Rules.vpfromvs(vs);
			den = Rules.denfromvp(vp);
			}
		/* Implement Q from Benites & Olsen, 2005 (WnLH & LA) */
		if (atten)
			{
			Qs = Rules.basinQs(vs);
			Qp = Rules.basinQp(Qs);
			}
		/* Implement minimum velocities */
		if (vp < Rules.minVp)
			vp = Rules.minVp;
		if (vs < Rules.minVs)
			vs = Rules.minVs;
		if (vs > Rules.maxv)
			vs = Rules.maxv;
		if (vp > Rules.maxv)
			vp = Rules.maxv;
		if (atten)
			return new ElastProps(vp, vs, den, Qp, Qs);
		return new ElastProps(vp, vs, den);
		}

	static ElastProps getRockGeotechProps(float dh, boolean atten, float z)
		{
		return getRockGeotechProps(zgeotech, rockgeotech, dh, atten, z);
		}
	static ElastProps getRockGeotechProps(float gtdepth, float rockvel, float dh, boolean atten, float z)
		{
		float den, vp, vs, vg, vb, Qs, Qp;
		Qs = Qp = -1F;
		if (z >= dh && z > gtdepth) /* geotech only affects surface of grid */
			{
			vp = Rules.bedvp(z);
			vs = Rules.vsfromvp(vp);
			den = Rules.denfromvp(vp);
			}
		else if (dh <= gtdepth) /* use only geotech velocity */
			{
			vs = rockvel;
			vp = Rules.vpfromvs(vs);
			den = Rules.denfromvp(vp);
			}
		else /* dh > gtdepth, use slowness-avg surface bedrock w/ geotech */
			{
			vg = gtdepth/rockvel;
			vb = (dh-gtdepth)/(Rules.vsfromvp(Rules.bedvp(0F)));
			vs = dh/(vg + vb);
			vp = Rules.vpfromvs(vs);
			den = Rules.denfromvp(vp);
			}
		/* Implement Q from Benites & Olsen, 2005 (WnLH & LA) */
		if (atten)
			{
			Qs = Rules.bedQs(vs);
			Qp = Rules.bedQp(Qs);
			}
		/* Implement minimum velocities */
		if (vp < Rules.minVp)
			vp = Rules.minVp;
		if (vs < Rules.minVs)
			vs = Rules.minVs;
		if (vs > Rules.maxv)
			vs = Rules.maxv;
		if (vp > Rules.maxv)
			vp = Rules.maxv;
		if (atten)
			return new ElastProps(vp, vs, den, Qp, Qs);
		return new ElastProps(vp, vs, den);
		}

	/* End of new Verison 5 Rules */

	/* Jachens and Blakely Great Basin sedimentary basin model */
	static float basinden(float z)
		{
		float den0=2.7F; /* g/cc */
		float dden[]={-0.75F, -0.65F, -0.55F, -0.35F, -0.25F}; /* g/cc */
		float dep[]={0F, 0.1F, 0.4F, 0.9F, 1.2F}; /* km */
		int nd=dden.length;
		int i;
		for (i=0; i<nd; i++)
			if (z <= dep[i])
				break;
		if (i==0)
			return den0 + dden[0];
		else if (i >= nd)
			return den0 + dden[nd-1];
		else
			return den0 + dden[i-1]
				+ (z - dep[i-1])*(dden[i] - dden[i-1])/(dep[i] - dep[i-1]);
		}
	/* Jachens and Blakely Great Basin volcanic rift and basin model */
	static float volcbasinden(float z)
		{
		float den0=2.7F; /* g/cc */
		float dden[]={-0.48F, -0.43F, -0.38F, -0.28F}; /* g/cc */
		float dep[]={0F, 0.2F, 0.6F, 1.2F}; /* km */
		int nd=dden.length;
		int i;
		for (i=0; i<nd; i++)
			if (z <= dep[i])
				break;
		if (i==0)
			return den0 + dden[0];
		else if (i >= nd)
			return den0 + dden[nd-1];
		else
			return den0 + dden[i-1]
				+ (z - dep[i-1])*(dden[i] - dden[i-1])/(dep[i] - dep[i-1]);
		}

	static boolean isValidGeotech(float vel) /* is a valid geotech Vs, in m/s */
		{
		return (vel > 0F && vel < 2000F);
		}

	static float getProximityThick(float dist)
		{
		float maxthick=6F; /* km */
		float maxdist=25F; /*km */
		if (dist <= 0F)
			return 0F;
		else if (dist >= maxdist)
			return maxthick;
		return maxthick*dist/maxdist;
		}

	static boolean isInBasin(Grid grid, float z, float thick)
		{
		/* Minimum basin thickness rule imposed only if thickness < dh */
		if (thick < grid.dh)
			return (z<thick && thick > Rules.minthick /*0.1*grid.dh*/);
		return z < thick;
		}
	/* Jachens geologic map unit values */
	static boolean isVolcanic(float geol)
		{
		return geol==400F;
		}
	static boolean isSedimentary(float geol)
		{
		return geol==300F;
		}
	static boolean isBedrock(float geol)
		{
		return geol==200F;
		}

	/* Amplifications computed from the geotech velocity and the basin depth at
	   each surface grid node, using the equation of E. H. Field (2000 SCEC Phase III
	   Report; 2001 USGS-OFR Poster). */
	static float getAmplif(float thick, float geotech)
		{
		return (float)(Math.pow(geotech, -0.704)*Math.exp(0.12*thick));
		}

	/***** Version 5 Additions to read a rules parameter line *****/
	static boolean readRulesLine(Vector input)
		{
		/* get grid parameter line and interpret */
		Rules.rulesline = new String(Rules.findRulesline(input));
		/* no rules parameter line, use NTS-LV defaults */
		if (Rules.rulesline == null || Rules.rulesline.length() < 5)
			return true;
		/*System.err.println(Rules.rulesline);*/
		if (!Rules.parseRulesLine())
			return false;
		System.err.println(Rules.makeInput());
		return true;
		}

	static boolean readRulesLine(String lines, CME cme)
		{
		/* get grid parameter line and interpret */
		Rules.rulesline = Rules.findRulesline(lines, cme);
		/* no rules parameter line, use NTS-LV defaults */
		if (Rules.rulesline == null || Rules.rulesline.length() < 5)
			return true;
		if (!Rules.parseRulesLine())
			return false;
		return true;
		}

	static String findRulesline(Vector input)
		{
		int nruleslines=0;
		String line, rulesline=null;
		for (int i=0; i<input.size(); i++)
			{
			line = (String)(input.elementAt(i));
			if (Getpar.isType(line, "rules"))
				{
				rulesline = new String(line);
				nruleslines++;
				}
			}
		/* Allow rules line to be left out, uses defaults
		if (nruleslines < 1)
			{
			System.err.println("ModelAssembler: Error: no rules definition lines in input, abort.");
			return null;
			}
		*/
		if (nruleslines > 1)
			{
			System.err.println("ModelAssembler: Warning: found >1="
			+ nruleslines + " rules-definition line(s) of infile; using last one in file.");
			}
		return rulesline;
		}

	static String findRulesline(String lines, CME cme)
		{
		String line, rulesline=null;
		rulesline = CME.getFirstLine(lines);
		System.err.println(rulesline);
		/* Allow rules line to be left out, uses defaults */
		if (!Getpar.isType(rulesline, "rules"))
			return null;
		return rulesline;
		}

	static String makeInput()
		{
		return "rules file=" + Rules.rulessource
		+ " Qf=" + PointVal.getSignif("" + Rules.Qf, 2)
		+ " minthick=" + PointVal.getSignif("" + Rules.minthick, 2)
		+ " vpovs=" + PointVal.getSignif("" + Rules.vpovs, 4)
		+ " minVs=" + PointVal.getSignif("" + Rules.minVs, 4)
		+ " minVp=" + PointVal.getSignif("" + Rules.minVp, 4)
		+ " maxv=" + PointVal.getSignif("" + Rules.maxv, 4)
		+ " zgeotech=" + PointVal.getSignif("" + Rules.zgeotech, 2)
		+ " rockgeotech=" + PointVal.getSignif("" + Rules.rockgeotech, 4)
		+ " soilgeotech=" + PointVal.getSignif("" + Rules.soilgeotech, 4);
		}

	static boolean parseRulesLine()
		{
		Getpar gp;
		gp = new Getpar(Rules.minthick);
		if (gp.getStringPar(Rules.rulesline, "file"))
			Rules.rulessource = gp.getString();
		if (gp.getFloatPar(Rules.rulesline, "minthick"))
			Rules.minthick = gp.getFloat();
		if (Rules.minthick <= 0F || Rules.minthick > 24000F)
			{
			System.err.println("ModelAssembler: Error: minthick="
			+ Rules.minthick + ", must be >0 and <24000 km, abort.");
			return false;
			}
		if (gp.getFloatPar(Rules.rulesline, "minVs"))
			Rules.minVs = gp.getFloat();
		if (Rules.minVs <= 0F || Rules.minVs > 20F)
			{
			System.err.println("ModelAssembler: Error: minVs="
			+ Rules.minVs + ", must be >0 and <20 km/s, abort.");
			return false;
			}
		if (gp.getFloatPar(Rules.rulesline, "minVp"))
			Rules.minVp = gp.getFloat();
		if (Rules.minVp <= 0F || Rules.minVp > 20F)
			{
			System.err.println("ModelAssembler: Error: minVp="
			+ Rules.minVp + ", must be >0 and <20 km/s, abort.");
			return false;
			}
		if (Rules.minVp < Rules.minVs)
			{
			System.err.println("ModelAssembler: Error: minVp="
			+ Rules.minVp + " is less than minVs="
			+ Rules.minVs + ", abort.");
			return false;
			}
		if (gp.getFloatPar(Rules.rulesline, "maxv"))
			Rules.maxv = gp.getFloat();
		if (Rules.maxv <= 0F || Rules.maxv > 20F)
			{
			System.err.println("ModelAssembler: Error: maxv="
			+ Rules.maxv + ", must be >0 and <20 km/s, abort.");
			return false;
			}
		if (gp.getFloatPar(Rules.rulesline, "vpovs"))
			Rules.vpovs = gp.getFloat();
		if (Rules.vpovs < 1F || Rules.vpovs > 1e20F)
			{
			System.err.println("ModelAssembler: Error: vpovs="
			+ Rules.vpovs + ", must be >=1 and <1e20, abort.");
			return false;
			}
		if (gp.getFloatPar(Rules.rulesline, "zgeotech"))
			Rules.zgeotech = gp.getFloat();
		if (Rules.zgeotech <= 0F || Rules.zgeotech > 24000F)
			{
			System.err.println("ModelAssembler: Error: zgeotech="
			+ Rules.zgeotech + ", must be >0 and <24000 km, abort.");
			return false;
			}
		if (gp.getFloatPar(Rules.rulesline, "rockgeotech"))
			Rules.rockgeotech = gp.getFloat();
		if (Rules.rockgeotech <= 0F || Rules.rockgeotech > 20F)
			{
			System.err.println("ModelAssembler: Error: rockgeotech="
			+ Rules.rockgeotech + ", must be >0 and <20 km/s, abort.");
			return false;
			}
		if (gp.getFloatPar(Rules.rulesline, "soilgeotech"))
			Rules.soilgeotech = gp.getFloat();
		if (Rules.soilgeotech <= 0F || Rules.soilgeotech > 20F)
			{
			System.err.println("ModelAssembler: Error: soilgeotech="
			+ Rules.soilgeotech + ", must be >0 and <20 km/s, abort.");
			return false;
			}
		
		return true;
		}

	} /* End of Rules class */