/* This file was generated by scm2java from source file "peaked-roof.scm" */
/*  "peaked-roof.scm" Calculate convection from peaked rectangular roofs. */
/*  Copyright (C) 2013 Aubrey Jaffer */

/*  This program is free software; you can redistribute it and/or modify */
/*  it under the terms of the GNU General Public License as published by */
/*  the Free Software Foundation, either version 3 of the License, or */
/*  (at your option) any later version. */
/*  */
/*  This program is distributed in the hope that it will be useful, */
/*  but WITHOUT ANY WARRANTY; without even the implied warranty of */
/*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the */
/*  GNU General Public License for more details. */
/*  */
/*  You should have received a copy of the GNU General Public License */
/*  along with this program; if not, write to the Free Software */
/*  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111, USA. */


import static java.lang.Math.*;

public class Peaked_roof {

    public static double lNorm(double x,double y,double p)
    {
	return pow((pow(x, p))+(pow(y, p)), 1/p);
    }

    public static double min(double x,double y,double z)
    {
	return Math.min(x, Math.min(y, z));
    }

    public static double max(double x,double y,double z)
    {
	return Math.max(x, Math.max(y, z));
    }

    public static double sq(double x)
    {
	return x*x;
    }


    public static double cbrt(double x)
    {
	return pow(x, 1/3.);
    }


    public static double nu_uscaa(double re,double pr,double l_eff,double eps__m)
    {
	return max((cbrt(pr))*(664.e-3*(sqrt(re))),
		   (cbrt(pr))*(37.e-3*(pow(re, 800.e-3))),
		   (re*(cbrt(pr)))/(sq(log10((eps__m/3.7)/l_eff)))/32);
    }

    /*  returns heat-flow-rate in W */
    /*  All angles are in degrees */

    public static double
	peakedRoofConvect(double h__m, // height
			  double w__m, // width
			  double eps__m, // mean height of surface roughness
			  double theta__deg, // inclination from vertical; negative is upward facing
			  double psi__deg, // forced flow angle from projection of vertical in plate
			  double v_f__m_per_s, // forced bulk flow velocity
			  double t_a__k, // ambient (bulk) air temperature
			  double t_s__k, // plate surface temperature
			  double p__pa,	 // bulk air pressure
			  double rh)	 // relative humidity
	throws Exception
    {

	double pi = 4 * atan(1);
	double pi_per_180 = pi / 180;
	double psi = pi_per_180 * psi__deg;
	double theta_e = pi_per_180 *
	    (t_a__k>t_s__k ? - (theta__deg) :theta__deg);
	double r_h = h__m / (abs(sin(theta_e)));
	double a = w__m * r_h;
	double q = w__m / r_h;
	double l = a / 2 / (r_h + w__m);
	double r = (Math.min(r_h, w__m)) / 2;
	double r_eff = r_h / 2;
	double g0 = 9.80665;
	double t__k = (t_s__k + t_a__k) / 2;
	double t__c = t__k - 273.15;
	double deltat__k = t_s__k - t_a__k;
	double t_a__c = t_a__k - 273.15;
	double p_v = exp(( -6.3536311e3 / t_a__k) +
			 34.04926034 +
			 (t_a__k * ( -19.509874e-3 + (t_a__k * 12.811805e-6))));
	double rho = ((p__pa - (rh * p_v)) / 287.058 / t__k) +
	    ((rh * p_v) / 461.495 / t__k);
	double m_a = 28.97;
	double m_v = 18.0153;
	double r_m = m_a / m_v;
	double c_pv = 1.869e3 + (t_a__c * ( -257.8e-3 + (t_a__c * 19.41e-3)));
	double c_pa = 1.034e3 + (t__c * ( -284.9e-3 + (t__c * 781.7e-6)));
	double x_v = (rh * p_v) / p__pa;
	double c_p = ((c_pa * (1 - x_v) * m_a) + (c_pv * x_v * m_v))  /
	    (((1 - x_v) * m_a) + (x_v * m_v));
	double x = (rh * p_v * m_v) / m_a / (p__pa - (rh * p_v));
	double mu_a1 = 16.74e-6 + ((t__c -  -10) * 49.25e-9);
	double mu_v1 = 9.217e-6 + (25.71e-9 * t__c);
	double mu =
	    (mu_a1 / (1 + (1.07 * x * r_m))) +
	    (mu_v1 / (1 + (930.e-3 / x / r_m)));
	double k_a1 = 22.41e-3 + (76.46e-6 * (t__k - 250));
	double k_v1 = 17.48e-3 + (69.46e-6 * t__c);
	double k =
	    (k_a1 / (1 + (1.07 * x * r_m))) +
	    (k_v1 / (1 + (930.e-3 / x / r_m)));
	double beta = 1. / t__k;
	double nu = mu / rho;
	double pr = (c_p * mu) / k;
	double ra =
	    (pr * (abs(sin(theta_e))) * beta * g0 *
	     (abs(deltat__k)) * (pow(l, 3))) /
	    (sq(nu));
	double ra_v =
	    (pr * (abs(cos(theta_e))) * beta * g0 *
	     (abs(deltat__k)) * (pow(r_eff, 3))) /
	    (sq(nu));
	double ra_r =
	    (pr * beta * g0 *
	     (abs(deltat__k)) * (pow(r, 3))) /
	    (sq(nu));
	double nu_9__3 =
	    sq(825.e-3 +
	       ((387.e-3 * (pow(ra_v, 166.66666666666666e-3))) /
		(pow(1 + (pow(492.e-3 / pr, 562.5e-3)),
		     296.2962962962963e-3))));
	double v_fe =
	    (nu * (min(sq(nu_9__3 / (cbrt(pr)) / 664.e-3),
		       pow(nu_9__3 / (cbrt(pr)) / 37.e-3, 1.25),
		       (nu_9__3 / (cbrt(pr))) *
		       (sq(log10((eps__m / 3.7) / r_eff))) * 32))) /
	    r_eff;
	double v_l =
	    sqrt((sq(v_f__m_per_s)) +
		 (2 * v_fe * v_f__m_per_s * (cos(psi))) +
		 (sq(v_fe)));
	double v_r =
	    sqrt((sq(v_f__m_per_s)) +
		 ( -2 * v_fe * v_f__m_per_s * (cos(psi))) +
		 (sq(v_fe)));
	double zeta =
	    atan2(v_f__m_per_s * (sin(psi)),
		  v_fe + (v_f__m_per_s * (cos(psi))));
	double l_eff =
	    (q>(abs(tan(zeta)))
	     ?(sqrt(a)) / (abs(cos(zeta))) / (sqrt(q)) /
	     (sq(1 + ((abs(tan(zeta))) / 3 / q)))
	     :(sqrt(a)) / (abs(sin(zeta))) / (pow(q,  -500.e-3)) /
	     (sq(1 + (q / 3 / (abs(tan(zeta)))))));
	double re_l = (v_l * l_eff) / nu;
	double re_r = (v_r * l_eff) / nu;
	if (eps__m>(100.e-3 * l_eff))
	    throw new
		Exception("peakedRoofConvect: eps__m too large: " + eps__m);
	return k * deltat__k * a *
	    (lNorm((0 < (sin(theta_e))
		    ?(544.e-3 * (pow(ra_r, 200.e-3))) /
		    (cbrt(1+(pow(785.e-3 / pr, 600.e-3)))) / r
		    :(sq(650.e-3 +
			 (360.e-3 * (pow(ra, 166.66666666666666e-3))))) / l),
		   ((nu_uscaa(re_l, pr, l_eff, eps__m)) +
		    (nu_uscaa(re_r, pr, l_eff, eps__m))) / 2 / l_eff,
		   4));
    }

    /*  returns h in W/(K.m^2) */
    /*  All angles are in degrees */

    public static double
	peakedRoofH(double h__m,
		    double w__m,
		    double eps__m,
		    double theta__deg,
		    double psi__deg,
		    double v_f__m_per_s,
		    double t_a__k,
		    double t_s__k,
		    double p__pa,
		    double rh)
	throws Exception
    {
	return (peakedRoofConvect(h__m,
				  w__m,
				  eps__m,
				  theta__deg,
				  psi__deg,
				  v_f__m_per_s,
				  t_a__k,
				  t_s__k,
				  p__pa,
				  rh)) /
	    h__m / w__m / (t_s__k - t_a__k);
    }

    /*  returns h in W/(m^2) */
    /*  All angles are in degrees */

    public static double
	peakedRoofHeatFlux(double h__m,
		    double w__m,
		    double eps__m,
		    double theta__deg,
		    double psi__deg,
		    double v_f__m_per_s,
		    double t_a__k,
		    double t_s__k,
		    double p__pa,
		    double rh)
	throws Exception
    {
	return (peakedRoofConvect(h__m,
				  w__m,
				  eps__m,
				  theta__deg,
				  psi__deg,
				  v_f__m_per_s,
				  t_a__k,
				  t_s__k,
				  p__pa,
				  rh)) /
	    h__m / w__m;
    }

    /*  returns W/K */
    /*  All angles are in degrees */

    public static double peakedRoofHa(double h__m,
				      double w__m,
				      double eps__m,
				      double theta__deg,
				      double psi__deg,
				      double v_f__m_per_s,
				      double t_a__k,
				      double t_s__k,
				      double p__pa,
				      double rh)
	throws Exception
    {
	return (peakedRoofConvect(h__m,
				  w__m,
				  eps__m,
				  theta__deg,
				  psi__deg,
				  v_f__m_per_s,
				  t_a__k,
				  t_s__k,
				  p__pa,
				  rh)) /
	    (t_s__k - t_a__k);
    }


    public static void main(String[] args) throws Exception
    {
	System.out.println((peakedRoofH(350.e-3, 350.e-3, 0,  -90, 0, 0, 294.15, 299.65, 101325, 500.e-3)));
	System.out.println((peakedRoofH(350.e-3, 350.e-3, 0,  -90, 0, 0, 299.65, 294.15, 101325, 500.e-3)));
	System.out.println();
	System.out.println((peakedRoofHa(300.e-3, 300.e-3, 3.e-3,  -90, 0, 1.02, 299.15, 294.15, 101325, 400.e-3)));
	System.out.println((peakedRoofHa(300.e-3, 300.e-3, 3.e-3,  -90, 0, 3.24, 299.15, 294.15, 101325, 400.e-3)));
	System.out.println((peakedRoofHa(300.e-3, 300.e-3, 3.e-3,  -90, 0, 10.2, 299.15, 294.15, 101325, 400.e-3)));
	System.out.println((peakedRoofHa(300.e-3, 300.e-3, 750.e-6,  -90, 0, 1.02, 299.15, 294.15, 101325, 400.e-3)));
	System.out.println((peakedRoofHa(300.e-3, 300.e-3, 750.e-6,  -90, 0, 3.24, 299.15, 294.15, 101325, 400.e-3)));
	System.out.println((peakedRoofHa(300.e-3, 300.e-3, 750.e-6,  -90, 0, 10.2, 299.15, 294.15, 101325, 400.e-3)));
	System.out.println((peakedRoofHa(300.e-3, 300.e-3, 0,  -90, 0, 1.02, 299.15, 294.15, 101325, 400.e-3)));
	System.out.println((peakedRoofHa(300.e-3, 300.e-3, 0,  -90, 0, 3.24, 299.15, 294.15, 101325, 400.e-3)));
	System.out.println((peakedRoofHa(300.e-3, 300.e-3, 0,  -90, 0, 10.2, 299.15, 294.15, 101325, 400.e-3)));
	return;
    }

    /*  These values match those calculated in section "Natural Convective */
    /*  Surface Conductance" to within 1%. */
    /*  */
    /*  4.059745840396925 */
    /*  1.250740104165605 */

    /*  These values match those calculated in section "Experiments" to */
    /*  within 1% */
    /*  */
    /*  1.7143255001959925	    5.445218742025886	    17.142101686932023 */
    /*  902.031242690098e-3	    2.864936666723055	    9.019105793138744 */
    /*  709.4462856503201e-3	    1.7880485697312907	    4.475240148242594 */


}