#!/usr/bin/perl

#
# Viestikallio 5m peilin tuulikuorman Az/El komponenttien
# estimointia eri elevaatiokulmilla ja tuulen azimuuteilla.
#
# Az = 0 on "suoraan edestä", az = 180 on "suoraan takaa".
# El = 0 on "peili suunnattu horisontiin", el = 90 on "suoraan ylös"
#

use Math::Trig;

# Taulukko joka kertoo kulman (asteina) vastaavuuden
# sivutuuliprofiilin poikkipinta-alaan prosentteina
@anglepercent = ();

# Työkalu jolla kyseinen epälineaarinen taulukko täytetään..
sub fill_anglepercent {
    local ($al0, $al1, $per0, $per1) = @_;
    local ($i, $j, $k, $x);

    $k = ($per1 - $per0) / ($al1 - $al0);

    for ($i = $al0; $i <= $al1; ++$i) {
	$j = 180 - $i;
	$x = $per0 + ($i - $al0) * $k;
	$anglepercent[$i] = $x;
	$anglepercent[$j] = $x;

	# printf "  %3d  %.3f\n", $i, $x;
    }
}

#%dishangles =
#    (  0 => 1.00,  15 => 0.97,
#       30 => 0.87, 45 => 0.72,
#       60 => 0.60, 65 => 0.58,
#       75 => 0.30, 80 => 0.25,
#       90 => 0.20 );

fill_anglepercent(   0,  15, 1.00, 0.97 );
fill_anglepercent(  15,  30, 0.97, 0.87 );
fill_anglepercent(  30,  45, 0.87, 0.72 );
fill_anglepercent(  45,  60, 0.72, 0.60 );
fill_anglepercent(  60,  75, 0.60, 0.30 );
fill_anglepercent(  75,  80, 0.30, 0.25 );
fill_anglepercent(  80,  90, 0.25, 0.20 );


$AZstep = 15;
$ELstep = 10;

#
# Ensin simppeleitä laskuja: otsapinnan suuruus..
#

printf "Peilin tuulta vasten olevan projektoidun otsapinnan suuruus %% (0/0 suunnassa: 20 m²)\n";
printf "AZ\\ EL";
for ($EL = 0; $EL < 91; $EL += $ELstep) {
    printf "  %5d", $EL;
}
printf "\n";
for ($AZ = 0; $AZ < 181; $AZ += $AZstep) {
    printf " %3d: ", $AZ;
    $az = deg2rad($AZ);

    for ($EL = 0; $EL < 91; $EL += $ELstep) {
	$el = deg2rad($EL);

	$a = cos($az)*cos($el);
	$p = int( rad2deg(acos(abs($a))) );
	printf "  %5.1f", $anglepercent[$p] * 100;

    } # EL

    printf "\n";
} # AZ
printf "\n";


printf "Peilin kokema tuulen suuntainen voima kN \@ 10 m/s\n";
printf "AZ\\ EL";
$F = 22.0/9.0;
for ($EL = 0; $EL < 91; $EL += $ELstep) {
    printf "  %5d", $EL;
}
printf "\n";
for ($AZ = 0; $AZ < 181; $AZ += $AZstep) {
    printf " %3d: ", $AZ;
    $az = deg2rad($AZ);

    for ($EL = 0; $EL < 91; $EL += $ELstep) {
	$el = deg2rad($EL);

	$a = cos($az)*cos($el);
	$p = int( rad2deg(acos(abs($a))) );
	printf "  %5.1f", $anglepercent[$p] * $F;

    } # EL

    printf "\n";
} # AZ
printf "\n";

printf "Peilin kokema tuulen suuntainen voima kN \@ 30 m/s\n";
printf "AZ\\ EL";
$F = 22.0/1.0;
for ($EL = 0; $EL < 91; $EL += $ELstep) {
    printf "  %5d", $EL;
}
printf "\n";
for ($AZ = 0; $AZ < 181; $AZ += $AZstep) {
    printf " %3d: ", $AZ;
    $az = deg2rad($AZ);

    for ($EL = 0; $EL < 91; $EL += $ELstep) {
	$el = deg2rad($EL);

	$a = cos($az)*cos($el);
	$p = int( rad2deg(acos(abs($a))) );
	printf "  %5.1f", $anglepercent[$p] * $F;

    } # EL

    printf "\n";
} # AZ
printf "\n";

printf "\n**** Elevaatioakselin laakeroinnin kokemat suorat voimat ****\n\n\n";

printf "Elevaatioakselin pitkittäinen tuulivoima kN \@ 10 m/s\n";
printf "AZ\\ EL";
$F = 22.0/9.0;
for ($EL = 0; $EL < 91; $EL += $ELstep) {
    printf "  %5d", $EL;
}
printf "\n";
for ($AZ = 0; $AZ < 181; $AZ += $AZstep) {
    printf " %3d: ", $AZ;
    $az = deg2rad($AZ);

    for ($EL = 0; $EL < 91; $EL += $ELstep) {
	$el = deg2rad($EL);

	$a = cos($az)*cos($el);
	$p = int( rad2deg(acos(abs($a))) );
	printf "  %5.1f", $anglepercent[$p] * sin($az) * $F;

    } # EL

    printf "\n";
} # AZ
printf "\n";

printf "Elevaatioakselin pitkittäinen tuulivoima kN \@ 30 m/s\n";
printf "AZ\\ EL";
$F = 22.0/1.0;
for ($EL = 0; $EL < 91; $EL += $ELstep) {
    printf "  %5d", $EL;
}
printf "\n";
for ($AZ = 0; $AZ < 181; $AZ += $AZstep) {
    printf " %3d: ", $AZ;
    $az = deg2rad($AZ);

    for ($EL = 0; $EL < 91; $EL += $ELstep) {
	$el = deg2rad($EL);

	$a = cos($az)*cos($el);
	$p = int( rad2deg(acos(abs($a))) );
	printf "  %5.1f", $anglepercent[$p] * sin($az) * $F;

    } # EL

    printf "\n";
} # AZ
printf "\n";



printf "Elevaatioakselin poikittainen tuulivoima kN \@ 10 m/s\n";
printf "AZ\\ EL";
$F = 22.0/9.0;
for ($EL = 0; $EL < 91; $EL += $ELstep) {
    printf "  %5d", $EL;
}
printf "\n";
for ($AZ = 0; $AZ < 181; $AZ += $AZstep) {
    printf " %3d: ", $AZ;
    $az = deg2rad($AZ);

    for ($EL = 0; $EL < 91; $EL += $ELstep) {
	$el = deg2rad($EL);

	$a = cos($az)*cos($el);
	$p = int( rad2deg(acos(abs($a))) );
	printf "  %5.1f", $anglepercent[$p] * cos($az) * $F;

    } # EL

    printf "\n";
} # AZ
printf "\n";

printf "Elevaatioakselin poikittainen tuulivoima kN \@ 30 m/s\n";
printf "AZ\\ EL";
$F = 22.0/1.0;
for ($EL = 0; $EL < 91; $EL += $ELstep) {
    printf "  %5d", $EL;
}
printf "\n";
for ($AZ = 0; $AZ < 181; $AZ += $AZstep) {
    printf " %3d: ", $AZ;
    $az = deg2rad($AZ);

    for ($EL = 0; $EL < 91; $EL += $ELstep) {
	$el = deg2rad($EL);

	$a = cos($az)*cos($el);
	$p = int( rad2deg(acos(abs($a))) );
	printf "  %5.1f", $anglepercent[$p] * cos($az) * $F;

    } # EL

    printf "\n";
} # AZ
printf "\n";


$X0 = 0.00; # Oikeakätisesti kohtisuorassa alempana tuleviin
$Y0 = 0.80; # Horisonttiin osoittavan peilin keilasuuntaan
$Z0 = 0.60; # Suoraan ylös


# Vaikeampia laskuja:  voimakeskiön asema akseliin nähden ja
# syntyvät momentit..


printf "\n\n***** Kiertomomentit azimuutin suhteen *****\n\n\n";

#$AZstep = 10;

printf "AZ akselin kiertomomentti kNm  \@ 10 m/s\n";
printf "AZ\\ EL";
$F = 22.0/9.0;
for ($EL = 0; $EL < 91; $EL += $ELstep) {
    printf "  %5d", $EL;
}
printf "\n";
for ($AZ = 0; $AZ < 181; $AZ += $AZstep) {
    printf " %3d: ", $AZ;
    $az = deg2rad($AZ);

    for ($EL = 0; $EL < 91; $EL += $ELstep) {
	$el = deg2rad($EL);

	$a = cos($az)*cos($el);
	$tuulikulma = int( rad2deg(acos(abs($a))) ); # Tuulen insidenssikulma
	$projektioala = $anglepercent[$tuulikulma];
	$tuulivoima = $projektioala * $F;

	# Lasketaan peilin voimakeskiön (X/Y/Z-0) sijainti suhteessa
	# tuulen suunnasta tulevaan maanpinnan suuntaiseen vektoriin

	# Kierretään peilin elevaatiota ylös, kaavassa kierretään
	# koordinatteja ALAS:
	$x1 = $X0;
	$y1 = $Y0 * cos(- $el) + $Z0 * sin(- $el);
	$z1 = $Z0 * cos(- $el) - $Y0 * sin(- $el);

	# printf "   %4.2f %4.2f %4.2f", $x1, $y1, $z1;

	# Kierretään peilin azimuuttia myötäpäivään, lasketaan
	# ikään kuin kierrettäisiin vastapäivään!
	# Kaavassa edelliseen nähden: z -> y, y -> x, (x -> z)

	$z2 = $z1;
	$x2 = $x1 * cos(- $az) + $y1 * sin(- $az);
	$y2 = $y1 * cos(- $az) - $x1 * sin(- $az);

	# printf "   %5.2f %5.2f %5.2f", $x2, $y2, $z2;

	# Nyt meillä on peilin keskiön kierretyt koordinaatit ja tuulivoima..

	# Az momentti on  X2 * tuulivoima !
	printf "  %5.1f", $x2 * $tuulivoima;

	# EL momentti on Z2 * tuulivoima * cos(AZ)
	#printf "  %5.1f", $z2 * cos($az) * $tuulivoima;

    } # EL

    printf "\n";
} # AZ
printf "\n";


printf "AZ akselin kiertomomentti kNm  \@ 30 m/s\n";
printf "AZ\\ EL";
$F = 22.0/1.0;
for ($EL = 0; $EL < 91; $EL += $ELstep) {
    printf "  %5d", $EL;
}
printf "\n";
for ($AZ = 0; $AZ < 181; $AZ += $AZstep) {
    printf " %3d: ", $AZ;
    $az = deg2rad($AZ);

    for ($EL = 0; $EL < 91; $EL += $ELstep) {
	$el = deg2rad($EL);

	$a = cos($az)*cos($el);
	$tuulikulma = int( rad2deg(acos(abs($a))) ); # Tuulen insidenssikulma
	$projektioala = $anglepercent[$tuulikulma];
	$tuulivoima = $projektioala * $F;

	# Lasketaan peilin voimakeskiön (X/Y/Z-0) sijainti suhteessa
	# tuulen suunnasta tulevaan maanpinnan suuntaiseen vektoriin

	# Kierretään peilin elevaatiota ylös, kaavassa kierretään
	# koordinatteja ALAS:
	$x1 = $X0;
	$y1 = $Y0 * cos(- $el) + $Z0 * sin(- $el);
	$z1 = $Z0 * cos(- $el) - $Y0 * sin(- $el);

	# printf "   %4.2f %4.2f %4.2f", $x1, $y1, $z1;

	# Kierretään peilin azimuuttia myötäpäivään, lasketaan
	# ikään kuin kierrettäisiin vastapäivään!
	# Kaavassa edelliseen nähden: z -> y, y -> x, (x -> z)

	$z2 = $z1;
	$x2 = $x1 * cos(- $az) + $y1 * sin(- $az);
	$y2 = $y1 * cos(- $az) - $x1 * sin(- $az);

	# printf "   %5.2f %5.2f %5.2f", $x2, $y2, $z2;

	# Nyt meillä on peilin keskiön kierretyt koordinaatit ja tuulivoima..

	# Az momentti on  X2 * tuulivoima !
	printf "  %5.1f", $x2 * $tuulivoima;

	# EL momentti on Z2 * tuulivoima * cos(AZ)
	#printf "  %5.1f", $z2 * cos($az) * $tuulivoima;

    } # EL

    printf "\n";
} # AZ
printf "\n";



printf "\n\n***** Kiertomomentit elevaation suhteen *****\n\n\n";


printf "EL akselin kiertomomentti kNm  \@ 10 m/s\n";
printf "AZ\\ EL";
$F = 22.0/9.0;
for ($EL = 0; $EL < 91; $EL += $ELstep) {
    printf "  %5d", $EL;
}
printf "\n";
for ($AZ = 0; $AZ < 181; $AZ += $AZstep) {
    printf " %3d: ", $AZ;
    $az = deg2rad($AZ);

    for ($EL = 0; $EL < 91; $EL += $ELstep) {
	$el = deg2rad($EL);

	$a = cos($az)*cos($el);
	$tuulikulma = int( rad2deg(acos(abs($a))) ); # Tuulen insidenssikulma
	$projektioala = $anglepercent[$tuulikulma];
	$tuulivoima = $projektioala * $F;

	# Lasketaan peilin voimakeskiön (X/Y/Z-0) sijainti suhteessa
	# tuulen suunnasta tulevaan maanpinnan suuntaiseen vektoriin

	# Kierretään peilin elevaatiota ylös, kaavassa kierretään
	# koordinatteja ALAS:
	$x1 = $X0;
	$y1 = $Y0 * cos(- $el) + $Z0 * sin(- $el);
	$z1 = $Z0 * cos(- $el) - $Y0 * sin(- $el);

	# printf "   %4.2f %4.2f %4.2f", $x1, $y1, $z1;

	# Kierretään peilin azimuuttia myötäpäivään, lasketaan
	# ikään kuin kierrettäisiin vastapäivään!
	# Kaavassa edelliseen nähden: z -> y, y -> x, (x -> z)

	$z2 = $z1;
	$x2 = $x1 * cos(- $az) + $y1 * sin(- $az);
	$y2 = $y1 * cos(- $az) - $x1 * sin(- $az);

	# printf "   %5.2f %5.2f %5.2f", $x2, $y2, $z2;

	# Nyt meillä on peilin keskiön kierretyt koordinaatit ja tuulivoima..

	# Az momentti on  X2 * tuulivoima !
	#printf "  %5.1f", $x2 * $tuulivoima;


	# EL momentti on Z2 * tuulivoima
	printf "  %5.1f", $z2 * cos($az) * $tuulivoima;

    } # EL

    printf "\n";
} # AZ
printf "\n";



printf "EL akselin kiertomomentti kNm  \@ 30 m/s\n";
printf "AZ\\ EL";
$F = 22.0/1.0;
for ($EL = 0; $EL < 91; $EL += $ELstep) {
    printf "  %5d", $EL;
}
printf "\n";
for ($AZ = 0; $AZ < 181; $AZ += $AZstep) {
    printf " %3d: ", $AZ;
    $az = deg2rad($AZ);

    for ($EL = 0; $EL < 91; $EL += $ELstep) {
	$el = deg2rad($EL);

	$a = cos($az)*cos($el);
	$tuulikulma = int( rad2deg(acos(abs($a))) ); # Tuulen insidenssikulma
	$projektioala = $anglepercent[$tuulikulma];
	$tuulivoima = $projektioala * $F;

	# Lasketaan peilin voimakeskiön (X/Y/Z-0) sijainti suhteessa
	# tuulen suunnasta tulevaan maanpinnan suuntaiseen vektoriin

	# Kierretään peilin elevaatiota ylös, kaavassa kierretään
	# koordinatteja ALAS:
	$x1 = $X0;
	$y1 = $Y0 * cos(- $el) + $Z0 * sin(- $el);
	$z1 = $Z0 * cos(- $el) - $Y0 * sin(- $el);

	# printf "   %4.2f %4.2f %4.2f", $x1, $y1, $z1;

	# Kierretään peilin azimuuttia myötäpäivään, lasketaan
	# ikään kuin kierrettäisiin vastapäivään!
	# Kaavassa edelliseen nähden: z -> y, y -> x, (x -> z)

	$z2 = $z1;
	$x2 = $x1 * cos(- $az) + $y1 * sin(- $az);
	$y2 = $y1 * cos(- $az) - $x1 * sin(- $az);

	# printf "   %5.2f %5.2f %5.2f", $x2, $y2, $z2;

	# Nyt meillä on peilin keskiön kierretyt koordinaatit ja tuulivoima..

	# Az momentti on  X2 * tuulivoima !
	#printf "  %5.1f", $x2 * $tuulivoima;


	# EL momentti on Z2 * tuulivoima
	printf "  %5.1f", $z2 * cos($az) * $tuulivoima;

    } # EL

    printf "\n";
} # AZ
printf "\n";