/* This file was generated by scm2java from source file "flat-roof.scm" */

// "Flat_roof.java" Calculate convection from flat 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 Flat_roof {

    public static double L_norm(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.);
    }

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

    public static double
	flatRoofConvect(double h__m, // height
			double w__m, // width
			double eps__m, // mean height of surface roughness
			double phi__deg, // rotation around center in plane of plate
			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 phi = pi_per_180 * phi__deg;
	double psi = pi_per_180 * psi__deg;
	double theta_e = pi_per_180 *
	    (t_a__k > t_s__k ? - theta__deg : theta__deg);
	double a = w__m * h__m;
	double q = w__m / h__m;
	double l = (h__m * w__m) / 2 / (h__m + w__m);
	double r = (Math.min(h__m, w__m)) / 2;
	double r_eff = (h__m * w__m) / (((abs(cos(phi))) * w__m) +
					((abs(sin(phi))) * h__m));
	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 =
	    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))))) -
	    phi;
	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 = (v * l_eff) / nu;
	double nu_uscaa =
	    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);

	if (eps__m > (100.e-3 * l_eff))
	    throw new Exception("flatRoofConvect: eps__m too large: " + eps__m);

	return
	    k * deltat__k * a *
	    L_norm((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 / l_eff,
		   4);
    }

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

    public static double flatRoofH(double h__m,
				   double w__m,
				   double eps__m,
				   double phi__deg,
				   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 (flatRoofConvect(h__m,
				w__m,
				eps__m,
				phi__deg,
				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 flatRoofHeatFlux(double h__m,
				   double w__m,
				   double eps__m,
				   double phi__deg,
				   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 (flatRoofConvect(h__m,
				w__m,
				eps__m,
				phi__deg,
				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 flatRoofHa(double h__m,
				    double w__m,
				    double eps__m,
				    double phi__deg,
				    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 (flatRoofConvect(h__m,
				w__m,
				eps__m,
				phi__deg,
				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((flatRoofH(350.e-3, 350.e-3, 0, 0,  -90, 0, 0, 294.15, 299.65, 101325, 500.e-3)));
	System.out.println((flatRoofH(350.e-3, 350.e-3, 0, 0,  -45, 0, 0, 294.15, 299.65, 101325, 500.e-3)));
	System.out.println((flatRoofH(350.e-3, 350.e-3, 0, 0, 0, 0, 0, 294.15, 299.65, 101325, 500.e-3)));
	System.out.println((flatRoofH(350.e-3, 350.e-3, 0, 0, 45, 0, 0, 294.15, 299.65, 101325, 500.e-3)));
	System.out.println((flatRoofH(350.e-3, 350.e-3, 0, 0, 90, 0, 0, 294.15, 299.65, 101325, 500.e-3)));
	System.out.println();
	System.out.println(flatRoofHa(300.e-3, 300.e-3, 3.e-3, 0, 90, 0, 1.02, 294.15, 299.15, 101325, 400.e-3));
	System.out.println(flatRoofHa(300.e-3, 300.e-3, 3.e-3, 0, 90, 0, 3.24, 294.15, 299.15, 101325, 400.e-3));
	System.out.println(flatRoofHa(300.e-3, 300.e-3, 3.e-3, 0, 90, 0, 10.2, 294.15, 299.15, 101325, 400.e-3));
	System.out.println(flatRoofHa(300.e-3, 300.e-3, 750.e-6, 0, 90, 0, 1.02, 294.15, 299.15, 101325, 400.e-3));
	System.out.println(flatRoofHa(300.e-3, 300.e-3, 750.e-6, 0, 90, 0, 3.24, 294.15, 299.15, 101325, 400.e-3));
	System.out.println(flatRoofHa(300.e-3, 300.e-3, 750.e-6, 0, 90, 0, 10.2, 294.15, 299.15, 101325, 400.e-3));
	System.out.println(flatRoofHa(300.e-3, 300.e-3, 0, 0, 90, 0, 1.02, 294.15, 299.15, 101325, 400.e-3));
	System.out.println(flatRoofHa(300.e-3, 300.e-3, 0, 0, 90, 0, 3.24, 294.15, 299.15, 101325, 400.e-3));
	System.out.println(flatRoofHa(300.e-3, 300.e-3, 0, 0, 90, 0, 10.2, 294.15, 299.15, 101325, 400.e-3));
	return;
    }

}

/* These values match those calculated in section "Natural Convective */
/* Surface Conductance" to within 1%. */
/* */
/* 4.059745840396925 */
/* 3.929299315446097 */
/* 2.9741925519602748 */
/* 2.721738472416408 */
/* 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 */