Monday, March 5, 2018

System File 5: A python script to calculate solar position



 A python script to calculate solar position without use of ephemeris library.
The script is imperfect. When the sun goes down, the position coordinates become imaginary numbers. Imaginary numbers are part of high level mathematics and not necessary for this project but are a default numerical type when certain conditions are present.
I may return to using an ephemeris library if the errors become unacceptable.

CODE
#NOAA Solar Calculations spreadsheet, Day

import math
import datetime
pi=3.14159
r=360/2/pi
B3=36.3 # Lat
B4=-82.2 # Lon
B5=TZ=-5 # time zone
#B not used
B7=11/26/2017
Y=datetime.datetime.now().year
M=datetime.datetime.now().month
D=datetime.datetime.now().day
H2="Blank"
# UTC-5=EDT
UTclock=datetime.datetime.utcnow().strftime("%H:%M:%S %Z")
LTclock=datetime.datetime.now().strftime("%H:%M:%S")
(hh, mm, ss)=UTclock.split(':')
UT=(int(hh)*3600+int(mm)*60+int(ss))/3600
LT=UT+TZ
h=datetime.datetime.now().hour
m=datetime.datetime.now().minute
s=datetime.datetime.now().second
if h > 12:
h=h-12
time=(int(h)*3600+int(m)*60+int(s))/3600
#Converting string into datetime
from datetime import datetime
datetime_object = datetime.strptime('Jun 1 2005  1:33PM', '%b %d %Y %I:%M%p')
# forward and reverse transformation
JD=(367 * Y ) - (int(7 * (Y + int((M + 9) / 12)) / 4)) + (int(275 * M / 9) + D + 1721013.5) + (UT / 24)
UTx=24 * ( JD + (int(7 * (Y + int((M + 9) / 12)) / 4)) - (int(275 * M / 9) + D + 1721013.5) - (367 * Y ) )
D2=int(JD)
E2=LT/24
F2=JD
#G=Julian Century
G2=(F2-2451545)/36525
#H=blank
#I=Geom Mean Long Sun#I2(1)=math.mod(280.46646+G2*(36000.76983+G2*0.0003032),360)
I2a=(280.46646+G2*(36000.76983+G2*0.0003032))
I2b=360
I2=I2a-int(I2a/I2b)*I2b
#J=Geometric Mean Anomaly of the Sun
J2=357.52911+G2*(35999.05029-0.0001537*G2)
#K=Eccentricity of Earth orbit
K2=0.016708634-G2*(0.000042037+0.0000001267*G2)
L2=math.sin(J2/r)*(1.914602-G2*(0.004817+0.000014*G2))+math.sin(2*J2/r)*(0.019993-0.000101*G2)+math.sin(3*J2/r)*0.000289
#M=Sun True Long degrees
M2=I2+L2
#N=Sun True Anomaly degrees
N2=J2+L2
#O=Sun Radius Vector, A.U.
O2=(1.000001018*(1-K2*K2))/(1+K2*math.cos(N2))
#P=Sun Apparent Long
P2=M2-0.00569-0.00478*math.sin(125.04-1934.136*G2)
#Q=Mean Oblique Ecliptic degrees (23.44)
Q2=23+(26+((21.448-G2*(46.815+G2*(0.00059-G2*0.001813))))/60)/60
#R=Oblique Correction factor (23.44)
R2=Q2+0.00256*math.cos((125.04-1934.136*G2)/r) # new
#S=Right Ascension degrees (120.74)
S2=math.atan2(math.cos(R2/r)*math.sin(P2/r), math.cos(P2/r))*r
T2=math.asin(math.sin(R2/r)*math.sin(P2/r))*r
#U=var y (0.04)
U2=math.tan(R2/2/r)*math.tan(R2/2/r)
#V=EoT minutes (13.55)
V2=4*(U2*math.sin(2*I2/r)-2*K2*math.sin(J2/r)+4*K2*U2*math.sin(J2/r)*math.cos(2*I2/r)-0.5*U2*U2*math.sin(4*I2/r)-1.25*K2*math.sin(2*J2/r))*r
V2=4*(U2*math.sin(2*(I2/r))-2*K2*math.sin((J2/r))+4*K2*U2*math.sin((J2/r))*math.cos(2*(I2/r))-0.5*U2*U2*math.sin(4*(I2/r))-1.25*K2*K2*math.sin(2*(J2/r)))*r
#W=HA Sunrise degrees (75.24)
W2=(math.acos(math.cos(90.833/r)/(math.cos(B3/r)*math.cos(T2/r))-math.tan(B3/r)*math.tan(T2/r)))*r#X=Sun noon LST (0.51)
X2=(720-4*B4-V2+B5*60)/1440
#Y=Sunrise time (0.3)
Y2=X2-W2*4/1440
#Z=Sunset time (0.72)
Z2=X2+W2*4/1440
#AA2=sunlight duration (601.9)
AA2=8*W2
#AB2=TST (404.75)
AB2a=(E2*1440+V2+4*B4-60*B5)
AB2b=1440
AB2=AB2a-int(AB2a/AB2b)*AB2b
if AB2/4 < 0 : AC2=AB2/4+180
else: AC2=AB2/4-180
#AD2=Zenith (93.46)
AD2=math.acos(math.sin(B3/r)*math.sin(T2/r)+math.cos(B3/r)*math.cos(T2/r)*math.cos(AC2/r))*r
#AE2= Sun elevation in degrees (-3.46)
AE2=90-AD2
#AF2=Atmospheric refraction in degrees (0.1)
if AE2>85 : AF2=0
elif AE2>5: AF2=58.1/math.tan((AE2))-0.07/math.pow(math.tan((AE2)),3)+0.000086/math.pow(math.tan((AE2)),5)
elif AE2>-0.575: AF2=1735+AE2*(-518.2+AE2*(103.4+AE2*(-12.79+AE2*0.711)))
else : AF2=-20.772/math.tan(AE2)/3600
 #AG2=Sun Elevation corrected for refraction (-3.37)
AG2=AE2+AF2
#AH2=Sun Azimuth degrees clockwise from North (112.93)
AH2a1=math.acos(((math.sin((B3/r)) * math.cos((AD2/r))) - math.sin((T2/r))) / ((math.cos((B3/r)) * math.sin(AD2/r))))*r + 180
AH2b1=360
AH2a2=540-math.acos(((math.sin(B3/r)*math.cos(AD2/r))-math.sin(T2/r))/((math.cos(B3/r)*math.sin(AD2/r))))*r
AH2b2=AH2b1
if AC2 > 0 : AH2=AH2a1-int(AH2a1/AH2b1)*AH2b1
else : AH2=AH2a2-int(AH2a2/AH2b2)*AH2b2
# error compensation
#AH2=AH2+7
'''
print(Y2+",")
print(Z2,",")
print(AE2,",")
print(AH2,",")
'''
#print (str(X2)+","+str(Y2)+","+str(AE2)+","+str(AH2));
ae2=str(AE2)
ah2=str(AH2)
#=============================================================================
# Si, Sm calculations
r=360/2/3.14159
a=0.14
h=0.535
ah=a*h
alpha=AE2 #40.5603 # Sel
beta=45
Z=90-alpha
psi=180
theta=AH2 #175.9508 # Saz

#AM=1/math.cos(Z/r) #  theta is 90-Sel or 90-alpha - use Z
#for flat earth, no earth curvature
AM=1/(math.cos(Z/r)+0.50572*pow(96.07995-Z/r,-1.6364))

# flat earth, no terrestrial elevation
#Si=1.353*math.pow(0.7,math.pow(AM,0.678))
Si=1.353*((1-ah)*pow(0.7,pow(AM,0.678))+ah)
Sm=Si*(math.cos(alpha/r)*math.sin(beta/r)*math.cos((psi-theta)/r)+math.sin(alpha/r)*math.cos(beta/r))
####################################################################################

print(str(ae2)+','+str(ah2)+','+str(Si)+','+str(Sm));
END CODE

System File 4: A javascript to handle date format differences between programming languages.

This java script makes sure the date is correctly interpreted between PHP and HTML. I found this on the web:

CODE
//<script>
/**
 * Return a formated string from a date Object mimicking PHP's date() functionality
 * format  string  "Y-m-d H:i:s" or similar PHP-style date format string
 * date    mixed   Date Object, Datestring, or milliseconds
 *
 * d-m-Y with leading zeroes
 * j-n-Y without leading zeroes
 *
 * dateFormat(d-m-Y, date), ex. date = 2-2-2017, and get 02-02-2017
 */
function dateFormat(format,date){
if(!date || date === "")
{ date = new Date();}
else if(typeof('date') !== 'object')
{date = new Date(date.replace(/-/g,"/")); // attempt to convert string to date object }
var string = '',
mo = date.getMonth(),   // month (0-11)
m1 = mo+1;     // month (1-12)
dow = date.getDay(),    // day of week (0-6)
d = date.getDate(),     // day of the month (1-31)
y = date.getFullYear(), // 1999 or 2003
h = date.getHours(),    // hour (0-23)
mi = date.getMinutes(), // minute (0-59)
s = date.getSeconds();  // seconds (0-59)
for (var i = 0, len = format.length; i < len; i++) {
switch(format[i])
{
case 'j': // Day of the month without leading zeros  (1 to 31)
string+= d;
break;
case 'd': // Day of the month, 2 digits WITH LEADING ZEROS (01 to 31)
string+= (d < 10) ? "0"+d : d;
break;
case 'm': // Numeric representation of a month, WITH LEADING ZEROS (01 to 12)
string+= (m1 < 10) ? "0"+m1 : m1;
break;
case 'n': // Numeric representation of a month, without leading zeros (1 to 12)
string+= m1;
break;
case 'Y': // A full numeric representation of a year, 4 digits (1999 OR 2003)
string+= y;
break;
default:
string+= format[i];
}
}
return string;
}
}
//</script>
END CODE

System File 3: HTML file to make an internet accessible real-time graph of panel activity


This file is a web page that graphs the current accumulated data for the day:

CODE
<!DOCTYPE html>
<html>
<head>
  <meta http-equiv="refresh" content="300" />
  <meta http-equiv="content-type" content="text/html; charset=UTF-8">
  <meta name="robots" content="noindex, nofollow">
  <meta name="googlebot" content="noindex, nofollow">
<script type="text/javascript" src="./resources/dateFormatscript.js"></script>
<script type="text/javascript" src="./resources/d3.v3.min.js"></script>
<script type="text/javascript" src="./resources/c3.js"></script>
<link rel="stylesheet" type="text/css" href="./resources/c3.css">
<style type="text/css"></style>
<title>once again</title>
<script type='text/javascript'>//<![CDATA[
window.onload=function(){
var d = new Date();
$d = (d.getMonth()+1) + '-' + d.getDate()+ '-' + d.getFullYear();
$de = 'data_' + dateFormat('m-d-Y',$d)+'.txt';
$df = dateFormat('m-d-Y',$d);

d3.csv($de)
.row(function(d) { return [d.time12, d.PVpower, d.BattSOC, d.elevcalc, d.azicalc]; })
.get(function(error,rows) {
rows.unshift (["time12","PVpower","BattSOC","elevcalc","azicalc"])
console.log(rows);
$c = new Date().getFullYear();
$a = new Date().getMonth() + 1;
$b = new Date().getUTCDate();
$d2 =  $df;
var chart = c3.generate({
title: {
text: $d2
},
  size:{
width: 1000,
height: 400
},
bindto: '#chart',
data: {
rows: rows,
type: 'line',
x: 'time12',
xFormat:'%I:%M:%S %p'
          },
point: {show: false},
tooltip: {show: false},
axis: {
x: {
localtime: true,
type:'timeseries',
tick:{
//culling:{max: 50},
fit: true,
count: 50,
format: '%I:%M:%S %p',
rotate: 65,
  }//end tick
},//end x

y: {
max: 350,
min: 0
} //end y
},//end axis
grid:{
x:{
show: true
},//end x
y:{
show: true
  }//end y
}//end grid
      });
});
}//]]>
</script>
</head>
<body>
  <div id="chart"></div>
<script>
  if (window.parent && window.parent.parent){
    window.parent.parent.postMessage(["resultsFrame", {
      height: document.body.getBoundingClientRect().height,
      slug: "None"
    }], "*")
  }
</script>
</body>
</html>
END CODE

System File 2:PHP script to write CSV header

File to write the CSV file header once each day or turn on cycle:

CODE
<?php
/*   epoutX-write-header.php   */

error_reporting(0);

exec ('sudo hwclock > /var/www/html/_hwc_out.txt');
exec ('date > /var/www/html/_date_out.txt');
exec ('sudo hwclock -w');

#$file= 'data_'.date('m-d-Y').'.txt'
$myfile = fopen('/var/www/html/data_'.date('m-d-Y').'.txt', "a") or die("Unable to open file!");
#shell_exec ( string 'sudo chown pi  $file'  )
fwrite($myfile,"time,time12,date,PV array voltage(V),PV array current(A),PVpower,Battery voltage(V),Battery charging current(A),Battery charging power(W),BattSOC,elevcalc,azicalc,Si,Sm\n");
fclose($myfile);

exec ('/bin/bash /var/www/html/epoutX-unlock_files.sh');



#$myfile = fopen('/var/www/html/data_XXX'.date('m-d-Y').'.txt', "a") or die("Unable to open file!");
#fwrite($myfile,"time,time12,date,PV array voltage(V),PV array current(A),PVpower,Battery voltage(V),Battery charging current(A),Battery charging power(W),BattSOC,azimuth,elevation\n");
#fclose($myfile);


?>
END CODE

System File 1:Data saver PHP script

 
   I'm publishing the files I've created to run and monitor, my system. There will be no explanations. I

The file to load and save data from the MPPT controller:

CODE
<?php
/*   epoutX-2.php   */
#date_default_timezone_set('EST');
#exec ('/bin/bash /var/www/html/epoutX-unlock_files.sh');

#$x=0;
require_once 'PhpEpsolarTracer.php';
$tracer = new PhpEpsolarTracer('/dev/ttyXRUSB0');
error_reporting(0);


$y = 1;
$x = 3;
#######################################
# exec ('/bin/bash /var/www/html/chd.sh');
#######################################

while($x = 3):
#######################################
exec ('/bin/bash /var/www/html/example_web.sh');
#######################################








#while($x=3):
$myfile = fopen('/var/www/html/data_'.date('m-d-Y').'.txt', "a") or die("Unable to open file!");
#If  ($y = 1){
# fwrite($myfile,"time,time12,date,PV array voltage(V),PV array current(A),PVpower,Battery voltage(V),Battery charging current(A),Battery charging power(W),BattSOC,azimuth,elevation,Si,Sm\n");
# $y = 0;
#  }


$nano = time_nanosleep(30, 0);
  if ($tracer->getRealtimeData()) {
    $i=0;
#--------------------------------------->
    $pythonnoaaephem = `/usr/bin/env python3.6 /var/www/html/resources/NOAA.py`;

#____________________________________


    $date=date("Y-m-d H:i:s");
    $seconds = strtotime($date);
    $seconds /= 30;
    $seconds = round($seconds);
    $seconds *= 30;
    $date2 = date("h:i:s", $seconds);
    $date = date("h:i:s A", $seconds);
    $today = date("Y-m-d"); # Y-m-d
    $now = $date;
#______________________________________
    fwrite($myfile, $date2);
    fwrite($myfile, ",");
    fwrite($myfile, $now);
    fwrite($myfile, ",");
    fwrite($myfile, $today.",");
#---------------------------------------->
    while ($i<=5):
      $dat[$i] = $tracer->realtimeData[$i].",";
       fwrite($myfile, $dat[$i]);
       $i++;
    endwhile;
  $i=12;
$dat[$i] = $tracer->realtimeData[$i].",";
    fwrite($myfile, $dat[$i]);

    }
#--------------------------------------------------------->
    fwrite($myfile,$pythonnoaaephem);
#--------------------------------------------------------->





endwhile;
    fclose($myfile);
?>
END CODE

Update and plans for CGI model of shadows



This post, I'm laying out plans to construct an animation of my property and the shadows cast throughout the year. I also recap progress on the installation.
   This will allow me to determine the best location for my array, or whether I will need to move the array, at least once per  year.

   The software I'll be using is Pov-Ray (Point of View Raytracer). If possible and necesasary, I'll use the architectural raytracer, Radiance.

   I have the formulae I need to do the calculations for solar positioning throughout the year. I'll get terrain data from D.E.M. data provided by the government for free. I also have twenty years of CDs containing GIS data and software to fall back on.
  The D.E.M. (Digital Elevation Map) is a file containing elevation data for points on an area of terrain. Generally, a seven degree quadrangle is the unit of measure. One 'quad' will contain data points similar to the pixels in an image file. The difference is the image data represents color or grayscale values. In fact, I have converted DEM files into grayscale high dynamic range images for use in rendering engines (raytracers). I constructed an animated gif of a local quad from topographic map graphic layers and a DEM layer. The gif shows the flat topo map flat, then the DEM data 'pushing' through from under to create a 3D map with a topo map overlay.

   I may use the Moray 3D object editing software to construct geometric tree and house objects to cast shadows. Some have created tree files for use in populating a landscape in Pov-Ray. I may use one of those.

   I've determined I need to double the number of watts for a comfortable minimum power for the system. A factor of two compensates for days that are overcast. There have to be more watts (panels) or more batteries. The panels last longer than batteries and are less expensive in terms of usage. Batteries last about five years. Panels have an industry wide warranty of twenty years for 80% rating.
   There is a balance between the number of batteries and panels.  Extra batteries are needed to draw only 50% of available power. The general rule for lead acid deep cycle batteries is: use no more than 50% to preserve longevity.

   There is one alternative to moving the array around. Find a way to get the owner of the adjoining property to remove or lower the top level of the evergreen trees blocking the sunlight from my array in the winter months.

   I may take that option. Roof mount would be an option if the roof were slanted more South. The roof points 20 degrees East of South. I would need to build a stage to position the aray pointing South. That is too difficult and costly. I don't like the idea of modifying my roof, either.

   That is the end of this post.