From 9dda5037d38af32f1b9a205d5227547ad69f1951 Mon Sep 17 00:00:00 2001 From: tmj-fstate Date: Tue, 28 Mar 2017 13:09:01 +0200 Subject: [PATCH] moon object, tweaks to the lighting model --- Globals.cpp | 5 +- Globals.h | 1 + Ground.cpp | 19 ++++- Model3d.cpp | 14 ++-- ResourceManager.cpp | 2 +- Segment.cpp | 4 +- World.cpp | 85 ++++++++++++++++----- World.h | 2 + maszyna.vcxproj.filters | 6 ++ skydome.cpp | 13 +++- sun.cpp | 161 +++++++++++++++++++++------------------- sun.h | 36 ++++----- 12 files changed, 220 insertions(+), 128 deletions(-) diff --git a/Globals.cpp b/Globals.cpp index f468c33d..5663318f 100644 --- a/Globals.cpp +++ b/Globals.cpp @@ -80,11 +80,12 @@ double Global::pCameraRotation; double Global::pCameraRotationDeg; std::vector Global::FreeCameraInit; std::vector Global::FreeCameraInitAngle; -double Global::fFogStart = 1700; -double Global::fFogEnd = 2000; float Global::Background[3] = {0.2f, 0.4f, 0.33f}; GLfloat Global::AtmoColor[] = {0.423f, 0.702f, 1.0f}; GLfloat Global::FogColor[] = {0.6f, 0.7f, 0.8f}; +double Global::fFogStart = 1700; +double Global::fFogEnd = 2000; +float Global::Overcast{ 0.1f }; // NOTE: all this weather stuff should be moved elsewhere #ifdef EU07_USE_OLD_LIGHTING_MODEL GLfloat Global::ambientDayLight[] = {0.40f, 0.40f, 0.45f, 1.0f}; // robocze GLfloat Global::diffuseDayLight[] = {0.55f, 0.54f, 0.50f, 1.0f}; diff --git a/Globals.h b/Globals.h index ba81d91e..26e9d8be 100644 --- a/Globals.h +++ b/Globals.h @@ -218,6 +218,7 @@ class Global static float Background[3]; static GLfloat AtmoColor[]; static GLfloat FogColor[]; + static float Overcast; // static bool bTimeChange; #ifdef EU07_USE_OLD_LIGHTING_MODEL static opengl_light AmbientLight; diff --git a/Ground.cpp b/Ground.cpp index df5cc8e6..7bfd2072 100644 --- a/Ground.cpp +++ b/Ground.cpp @@ -2798,16 +2798,29 @@ bool TGround::Init(std::string File) { // Ra: ustawienie parametrów OpenGL przeniesione do FirstInit WriteLog("Scenery atmo definition"); parser.getTokens(3); - parser >> Global::AtmoColor[0] >> Global::AtmoColor[1] >> Global::AtmoColor[2]; + parser + >> Global::AtmoColor[0] + >> Global::AtmoColor[1] + >> Global::AtmoColor[2]; parser.getTokens(2); - parser >> Global::fFogStart >> Global::fFogEnd; + parser + >> Global::fFogStart + >> Global::fFogEnd; if (Global::fFogEnd > 0.0) { // ostatnie 3 parametry są opcjonalne parser.getTokens(3); - parser >> Global::FogColor[0] >> Global::FogColor[1] >> Global::FogColor[2]; + parser + >> Global::FogColor[0] + >> Global::FogColor[1] + >> Global::FogColor[2]; } parser.getTokens(); parser >> token; + if( token != "endatmo" ) { + // optional overcast parameter + // NOTE: parameter system needs some decent replacement, but not worth the effort if we're moving to built-in editor + Global::Overcast = clamp( std::stof( token ), 0.0f, 1.0f ); + } while (token.compare("endatmo") != 0) { // a kolejne parametry są pomijane parser.getTokens(); diff --git a/Model3d.cpp b/Model3d.cpp index 4d9fb42d..efd08a4b 100644 --- a/Model3d.cpp +++ b/Model3d.cpp @@ -313,12 +313,16 @@ int TSubModel::Load(cParser &parser, TModel3d *Model, int Pos, bool dynamic) } std::string discard; parser.getTokens(13, false); - parser >> fNearAttenStart >> discard >> fNearAttenEnd >> discard >> bUseNearAtten >> - discard >> iFarAttenDecay >> discard >> fFarDecayRadius >> discard >> - fCosFalloffAngle // kąt liczony dla średnicy, a nie promienia + parser + >> fNearAttenStart + >> discard >> fNearAttenEnd + >> discard >> bUseNearAtten + >> discard >> iFarAttenDecay + >> discard >> fFarDecayRadius + >> discard >> fCosFalloffAngle // kąt liczony dla średnicy, a nie promienia >> discard >> fCosHotspotAngle; // kąt liczony dla średnicy, a nie promienia - fCosFalloffAngle = cos(DegToRad(0.5 * fCosFalloffAngle)); - fCosHotspotAngle = cos(DegToRad(0.5 * fCosHotspotAngle)); + fCosFalloffAngle = std::cos( DegToRad( 0.5f * fCosFalloffAngle ) ); + fCosHotspotAngle = std::cos( DegToRad( 0.5f * fCosHotspotAngle ) ); iNumVerts = 1; /* iFlags |= 0x4010; // rysowane w cyklu nieprzezroczystych, macierz musi zostać bez zmiany diff --git a/ResourceManager.cpp b/ResourceManager.cpp index 412aa853..12d44b38 100644 --- a/ResourceManager.cpp +++ b/ResourceManager.cpp @@ -18,7 +18,7 @@ double ResourceManager::_lastReport = 0.0f; void ResourceManager::Register(Resource *resource) { - _resources.push_back(resource); + _resources.emplace_back(resource); }; void ResourceManager::Unregister(Resource *resource) diff --git a/Segment.cpp b/Segment.cpp index 821bb407..3106ee28 100644 --- a/Segment.cpp +++ b/Segment.cpp @@ -401,7 +401,9 @@ void TSegment::RenderLoft(const vector6 *ShapePoints, int iNumShapePoints, doubl jmm1 * ShapePoints[j].z + m1 * ShapePoints[j + iNumShapePoints].z, tv1); glVertex3f(pt.x, pt.y, pt.z); // pt nie mamy gdzie zapamiętać? } - if (p) // jeśli jest wskaźnik do tablicy + // BUG: things blow up badly in the following part in 64bit version on baltyk.scn + // TODO: sort this mess out when the time comes to reorganize spline generation + if( p ) // jeśli jest wskaźnik do tablicy if (*p) if (!j) // to dla pierwszego punktu { diff --git a/World.cpp b/World.cpp index 3b1d837d..febe4749 100644 --- a/World.cpp +++ b/World.cpp @@ -2311,37 +2311,75 @@ void world_environment::init() { m_sun.init(); + m_moon.init(); m_stars.init(); m_clouds.Init(); } void world_environment::update() { - - // move sun... + // move celestial bodies... m_sun.update(); - auto const position = m_sun.getPosition(); - // ...update the global data to match new sun state... - Global::SunAngle = m_sun.getAngle(); - Global::DayLight.set_position( position ); - Global::DayLight.direction = -1.0 * m_sun.getDirection(); + m_moon.update(); + // ...determine source of key light and adjust global state accordingly... + auto const sunlightlevel = m_sun.getIntensity(); + auto const moonlightlevel = m_moon.getIntensity(); + float keylightintensity; + float twilightfactor; + float3 keylightcolor; + if( moonlightlevel > sunlightlevel ) { + // rare situations when the moon is brighter than the sun, typically at night + Global::SunAngle = m_moon.getAngle(); + Global::DayLight.set_position( m_moon.getPosition() ); + Global::DayLight.direction = -1.0 * m_moon.getDirection(); + keylightintensity = moonlightlevel; + // if the moon is up, it overrides the twilight + twilightfactor = 0.0f; + keylightcolor = float3( 255.0f / 255.0f, 242.0f / 255.0f, 202.0f / 255.0f ); + } + else { + // regular situation with sun as the key light + Global::SunAngle = m_sun.getAngle(); + Global::DayLight.set_position( m_sun.getPosition() ); + Global::DayLight.direction = -1.0 * m_sun.getDirection(); + keylightintensity = sunlightlevel; + // diffuse (sun) intensity goes down after twilight, and reaches minimum 18 degrees below horizon + twilightfactor = clamp( -Global::SunAngle, 0.0f, 18.0f ) / 18.0f; + // TODO: crank orange up at dawn/dusk + keylightcolor = float3( 255.0f / 255.0f, 242.0f / 255.0f, 231.0f / 255.0f ); + float const duskfactor = 1.0f - clamp( Global::SunAngle, 0.0f, 18.0f ) / 18.0f; + keylightcolor = interpolate( + float3( 255.0f / 255.0f, 242.0f / 255.0f, 231.0f / 255.0f ), + float3( 235.0f / 255.0f, 140.0f / 255.0f, 36.0f / 255.0f ), + duskfactor ); + } // ...update skydome to match the current sun position as well... - m_skydome.Update( position ); + m_skydome.SetOvercastFactor( Global::Overcast ); + m_skydome.Update( m_sun.getPosition() ); // ...retrieve current sky colour and brightness... auto const skydomecolour = m_skydome.GetAverageColor(); auto const skydomehsv = RGBtoHSV( skydomecolour ); - auto const intensity = std::min( 1.15f * (0.05f + m_sun.getIntensity() + skydomehsv.z), 1.25f ); - // ...update light colours and intensity. - // NOTE: intensity combines intensity of the sun and the light reflected by the sky dome + // sun strength is reduced by overcast level + keylightintensity *= ( 1.0f - Global::Overcast * 0.65f ); + + // intensity combines intensity of the sun and the light reflected by the sky dome // it'd be more technically correct to have just the intensity of the sun here, // but whether it'd _look_ better is something to be tested - Global::DayLight.diffuse[ 0 ] = intensity * 255.0f / 255.0f; - Global::DayLight.diffuse[ 1 ] = intensity * 242.0f / 255.0f; - Global::DayLight.diffuse[ 2 ] = intensity * 231.0f / 255.0f; + auto const intensity = std::min( 1.15f * ( 0.05f + keylightintensity + skydomehsv.z ), 1.25f ); + // the impact of sun component is reduced proportionally to overcast level, as overcast increases role of ambient light + auto const diffuselevel = interpolate( keylightintensity, intensity * ( 1.0f - twilightfactor ), 1.0f - Global::Overcast * 0.75f ); + // ...update light colours and intensity. + keylightcolor = keylightcolor * diffuselevel; + Global::DayLight.diffuse[ 0 ] = keylightcolor.x; + Global::DayLight.diffuse[ 1 ] = keylightcolor.y; + Global::DayLight.diffuse[ 2 ] = keylightcolor.z; - Global::DayLight.ambient[ 0 ] = skydomecolour.x; - Global::DayLight.ambient[ 1 ] = skydomecolour.y; - Global::DayLight.ambient[ 2 ] = skydomecolour.z; + // tonal impact of skydome color is inversely proportional to how high the sun is above the horizon + // (this is pure conjecture, aimed more to 'look right' than be accurate) + float const ambienttone = clamp( 1.0f - ( Global::SunAngle / 90.0f ), 0.0f, 1.0f ); + Global::DayLight.ambient[ 0 ] = interpolate( skydomehsv.z, skydomecolour.x, ambienttone ); + Global::DayLight.ambient[ 1 ] = interpolate( skydomehsv.z, skydomecolour.y, ambienttone ); + Global::DayLight.ambient[ 2 ] = interpolate( skydomehsv.z, skydomecolour.z, ambienttone ); Global::fLuminance = intensity; @@ -2379,11 +2417,22 @@ world_environment::render() { m_skydome.Render(); m_stars.render(); - m_clouds.Render( m_skydome.GetAverageColor() * 2.5f ); + float const duskfactor = 1.0f - clamp( std::abs(m_sun.getAngle()), 0.0f, 12.0f ) / 12.0f; + float3 suncolor = interpolate( + m_skydome.GetAverageColor(), + float3( 235.0f / 255.0f, 140.0f / 255.0f, 36.0f / 255.0f ), + duskfactor ) * m_skydome.GetAverageColor().z; + + m_clouds.Render( + interpolate( m_skydome.GetAverageColor(), suncolor, duskfactor * 0.65f ) +// m_skydome.GetAverageColor() + * ( 1.0f - Global::Overcast * 0.5f ) // overcast darkens the clouds + * 2.5f ); // arbitrary adjustment factor if( DebugModeFlag == true ) { // mark sun position for easier debugging m_sun.render(); + m_moon.render(); } Global::DayLight.apply_angle(); Global::DayLight.apply_intensity(); diff --git a/World.h b/World.h index 551f8fed..119f62ad 100644 --- a/World.h +++ b/World.h @@ -15,6 +15,7 @@ http://mozilla.org/MPL/2.0/. #include "Ground.h" #include "sky.h" #include "sun.h" +#include "moon.h" #include "stars.h" #include "skydome.h" #include "mczapkie/mover.h" @@ -33,6 +34,7 @@ private: CSkyDome m_skydome; cStars m_stars; cSun m_sun; + cMoon m_moon; TSky m_clouds; }; diff --git a/maszyna.vcxproj.filters b/maszyna.vcxproj.filters index ba97afe4..9ccc0f0b 100644 --- a/maszyna.vcxproj.filters +++ b/maszyna.vcxproj.filters @@ -210,6 +210,9 @@ Source Files + + Source Files + @@ -407,6 +410,9 @@ Header Files + + Header Files + diff --git a/skydome.cpp b/skydome.cpp index 2afb16c3..27074616 100644 --- a/skydome.cpp +++ b/skydome.cpp @@ -184,7 +184,7 @@ void CSkyDome::SetGammaCorrection( float const Gamma ) { void CSkyDome::SetOvercastFactor( float const Overcast ) { - m_overcast = clamp( Overcast, 0.0f, 1.0f ); + m_overcast = clamp( Overcast, 0.0f, 1.0f ) * 0.75f; // going above 0.65 makes the model go pretty bad, appearance-wise } void CSkyDome::GetPerez( float *Perez, float Distribution[ 5 ][ 2 ], const float Turbidity ) { @@ -277,7 +277,7 @@ void CSkyDome::RebuildColors() { // luminance(Y) for clear & overcast sky float const yclear = PerezFunctionO2( perezluminance, icostheta, gamma, cosgamma2, zenithluminance ); - float const yover = ( 1.0f + 2.0f * vertex.y ) / 3.0f; + float const yover = zenithluminance * ( 1.0f + 2.0f * vertex.y ) / 3.0f; float const Y = interpolate( yclear, yover, m_overcast ); float const X = (x / y) * Y; @@ -295,6 +295,11 @@ void CSkyDome::RebuildColors() { colorconverter.z = 1.0f - exp( -m_expfactor * colorconverter.z ); } + // desaturate sky colour, based on overcast level + if( colorconverter.y > 0.0f ) { + colorconverter.y *= ( 1.0f - m_overcast ); + } + // override the hue, based on sun height above the horizon. crude way to deal with model shortcomings // correction begins when the sun is higher than 10 degrees above the horizon, and fully in effect at 10+15 degrees float const degreesabovehorizon = 90.0f - m_thetasun * ( 180.0f / M_PI ); @@ -318,8 +323,8 @@ void CSkyDome::RebuildColors() { */ // crude correction for the times where the model breaks (late night) // TODO: use proper night sky calculation for these times instead - if( ( color.x <= 0.0f ) - && ( color.y <= 0.0f ) ) { + if( ( color.x <= 0.05f ) + && ( color.y <= 0.05f ) ) { // darken the sky as the sun goes deeper below the horizon // 15:50:75 is picture-based night sky colour. it may not be accurate but looks 'right enough' color.z = 0.75f * std::max( color.z + m_sundirection.y, 0.075f ); diff --git a/sun.cpp b/sun.cpp index 1d0fbfea..6d671212 100644 --- a/sun.cpp +++ b/sun.cpp @@ -48,8 +48,6 @@ cSun::render() { GLfloat LightPosition[]= { 10.0f, 50.0f, -5.0f, 1.0f }; // ambient glLightfv(GL_LIGHT1, GL_POSITION, LightPosition ); */ - glDisable(GL_LIGHTING); - glDisable(GL_FOG); glColor4f( 255.0f/255.0f, 242.0f/255.0f, 231.0f/255.0f, 1.f ); // debug line to locate the sun easier Math3D::vector3 position = m_position; @@ -62,8 +60,6 @@ cSun::render() { // radius is a result of scaling true distance down to 2km -- it's scaled by equal ratio gluSphere( sunsphere, (float)(m_body.distance * 9.359157), 12, 12 ); glPopMatrix(); - glEnable(GL_FOG); - glEnable(GL_LIGHTING); } Math3D::vector3 @@ -115,106 +111,117 @@ void cSun::setPressure( float const Pressure ) { void cSun::move() { - static double degrad = 57.295779513; // converts from radians to degrees - static double raddeg = 0.0174532925; // converts from degrees to radians + static double radtodeg = 57.295779513; // converts from radians to degrees + static double degtorad = 0.0174532925; // converts from degrees to radians SYSTEMTIME localtime; // time for the calculation time( &localtime ); - if( m_observer.hour >= 0 ) { localtime.wHour = m_observer.hour; } + if( m_observer.hour >= 0 ) { localtime.wHour = m_observer.hour; } if( m_observer.minute >= 0 ) { localtime.wMinute = m_observer.minute; } if( m_observer.second >= 0 ) { localtime.wSecond = m_observer.second; } double ut = localtime.wHour - + localtime.wMinute / 60.0 // too low resolution, noticeable skips - + localtime.wSecond / 3600.0; // good enough in normal circumstances -/* - + localtime.wMilliseconds / 3600000.0; // for really smooth movement -*/ - double daynumber = 367 * localtime.wYear - - 7 * ( localtime.wYear + ( localtime.wMonth + 9 ) /12 ) / 4 - + 275 * localtime.wMonth / 9 - + localtime.wDay - - 730530 - + (ut / 24.0); + + localtime.wMinute / 60.0 // too low resolution, noticeable skips + + localtime.wSecond / 3600.0; // good enough in normal circumstances + /* + + localtime.wMilliseconds / 3600000.0; // for really smooth movement + */ + double daynumber = 367 * localtime.wYear + - 7 * ( localtime.wYear + ( localtime.wMonth + 9 ) / 12 ) / 4 + + 275 * localtime.wMonth / 9 + + localtime.wDay + - 730530 + + ( ut / 24.0 ); - // Universal Coordinated (Greenwich standard) time + // Universal Coordinated (Greenwich standard) time m_observer.utime = ut * 3600.0; m_observer.utime = m_observer.utime / 3600.0 - m_observer.timezone; - - // mean longitude - m_body.mnlong = 280.460 + 0.9856474 * daynumber; - m_body.mnlong -= 360.0 * (int) ( m_body.mnlong / 360.0 ); // clamp the range to 0-360 - if( m_body.mnlong < 0.0 ) m_body.mnlong += 360.0; - - // mean anomaly - m_body.mnanom = 357.528 + 0.9856003 * daynumber; - m_body.mnanom -= 360.0 * (int) ( m_body.mnanom / 360.0 ); // clamp the range to 0-360 - if( m_body.mnanom < 0.0 ) m_body.mnanom += 360.0; - - // ecliptic longitude - m_body.eclong = m_body.mnlong - + 1.915 * sin( m_body.mnanom * raddeg ) - + 0.020 * sin ( 2.0 * m_body.mnanom * raddeg ); - m_body.eclong -= 360.0 * (int)( m_body.eclong / 360.0 ); - if( m_body.eclong < 0.0 ) m_body.eclong += 360.0; // clamp the range to 0-360 - + // perihelion longitude + m_body.phlong = 282.9404 + 4.70935e-5 * daynumber; // w + // orbit eccentricity + double const e = 0.016709 - 1.151e-9 * daynumber; + // mean anomaly + m_body.mnanom = 356.0470 + 0.9856002585 * daynumber; // M + m_body.mnanom -= 360.0 * (int)( m_body.mnanom / 360.0 ); // clamp the range to 0-360 + if( m_body.mnanom < 0.0 ) m_body.mnanom += 360.0; // obliquity of the ecliptic - m_body.ecobli = 23.439 - 4.0e-07 * daynumber; + m_body.oblecl = 23.4393 - 3.563e-7 * daynumber; + // mean longitude + m_body.mnlong = m_body.phlong + m_body.mnanom; // L = w + M + m_body.mnlong -= 360.0 * (int)( m_body.mnlong / 360.0 ); // clamp the range to 0-360 + if( m_body.mnlong < 0.0 ) m_body.mnlong += 360.0; + // eccentric anomaly + double const E = m_body.mnanom + radtodeg * e * std::sin( degtorad * m_body.mnanom ) * ( 1.0 + e * std::cos( degtorad * m_body.mnanom ) ); + // ecliptic plane rectangular coordinates + double const xv = std::cos( degtorad * E ) - e; + double const yv = std::sin( degtorad * E ) * std::sqrt( 1.0 - e*e ); + // distance + m_body.distance = std::sqrt( xv*xv + yv*yv ); // r + // true anomaly + m_body.tranom = radtodeg * std::atan2( yv, xv ); // v + // ecliptic longitude + m_body.eclong = m_body.tranom + m_body.phlong; // lon = v + w + m_body.eclong -= 360.0 * (int)( m_body.eclong / 360.0 ); + if( m_body.eclong < 0.0 ) m_body.eclong += 360.0; // clamp the range to 0-360 +/* + // ecliptic rectangular coordinates + double const x = m_body.distance * std::cos( degtorad * m_body.eclong ); + double const y = m_body.distance * std::sin( degtorad * m_body.eclong ); + double const z = 0.0; + // equatorial rectangular coordinates + double const xequat = x; + double const yequat = y * std::cos( degtorad * m_body.oblecl ) - 0.0 * std::sin( degtorad * m_body.oblecl ); + double const zequat = y * std::sin( degtorad * m_body.oblecl ) + 0.0 * std::cos( degtorad * m_body.oblecl ); + // declination + m_body.declin = radtodeg * std::atan2( zequat, std::sqrt( xequat*xequat + yequat*yequat ) ); + // right ascension + m_body.rascen = radtodeg * std::atan2( yequat, xequat ); +*/ + // declination + m_body.declin = radtodeg * std::asin( std::sin( m_body.oblecl * degtorad ) * std::sin( m_body.eclong * degtorad ) ); - // declination - m_body.declin = degrad * asin( sin (m_body.ecobli * raddeg) * sin (m_body.eclong * raddeg) ); + // right ascension + double top = std::cos( degtorad * m_body.oblecl ) * std::sin( degtorad * m_body.eclong ); + double bottom = std::cos( degtorad * m_body.eclong ); - // right ascension - double top = cos ( raddeg * m_body.ecobli ) * sin ( raddeg * m_body.eclong ); - double bottom = cos ( raddeg * m_body.eclong ); + m_body.rascen = radtodeg * std::atan2( top, bottom ); + if( m_body.rascen < 0.0 ) m_body.rascen += 360.0; // (make it a positive angle) - m_body.rascen = degrad * atan2( top, bottom ); - if( m_body.rascen < 0.0 ) m_body.rascen += 360.0; // (make it a positive angle) + // Greenwich mean sidereal time + m_observer.gmst = 6.697375 + 0.0657098242 * daynumber + m_observer.utime; - // Greenwich mean sidereal time - m_observer.gmst = 6.697375 + 0.0657098242 * daynumber + m_observer.utime; + m_observer.gmst -= 24.0 * (int)( m_observer.gmst / 24.0 ); + if( m_observer.gmst < 0.0 ) m_observer.gmst += 24.0; - m_observer.gmst -= 24.0 * (int)( m_observer.gmst / 24.0 ); - if( m_observer.gmst < 0.0 ) m_observer.gmst += 24.0; + // local mean sidereal time + m_observer.lmst = m_observer.gmst * 15.0 + m_observer.longitude; - // local mean sidereal time - m_observer.lmst = m_observer.gmst * 15.0 + m_observer.longitude; + m_observer.lmst -= 360.0 * (int)( m_observer.lmst / 360.0 ); + if( m_observer.lmst < 0.0 ) m_observer.lmst += 360.0; - m_observer.lmst -= 360.0 * (int)( m_observer.lmst / 360.0 ); - if( m_observer.lmst < 0.0 ) m_observer.lmst += 360.0; + // hour angle + m_body.hrang = m_observer.lmst - m_body.rascen; - // hour angle - m_body.hrang = m_observer.lmst - m_body.rascen; - - if( m_body.hrang < -180.0 ) m_body.hrang += 360.0; // (force it between -180 and 180 degrees) + if( m_body.hrang < -180.0 ) m_body.hrang += 360.0; // (force it between -180 and 180 degrees) else if( m_body.hrang > 180.0 ) m_body.hrang -= 360.0; - double cz; // cosine of the solar zenith angle + double cz; // cosine of the solar zenith angle - double tdatcd = cos( raddeg * m_body.declin ); - double tdatch = cos( raddeg * m_body.hrang ); - double tdatcl = cos( raddeg * m_observer.latitude ); - double tdatsd = sin( raddeg * m_body.declin ); - double tdatsl = sin( raddeg * m_observer.latitude ); + double tdatcd = std::cos( degtorad * m_body.declin ); + double tdatch = std::cos( degtorad * m_body.hrang ); + double tdatcl = std::cos( degtorad * m_observer.latitude ); + double tdatsd = std::sin( degtorad * m_body.declin ); + double tdatsl = std::sin( degtorad * m_observer.latitude ); - cz = tdatsd * tdatsl + tdatcd * tdatcl * tdatch; + cz = tdatsd * tdatsl + tdatcd * tdatcl * tdatch; // (watch out for the roundoff errors) - if( fabs (cz) > 1.0 ) { cz >= 0.0 ? cz = 1.0 : cz = -1.0; } + if( fabs( cz ) > 1.0 ) { cz >= 0.0 ? cz = 1.0 : cz = -1.0; } - m_body.zenetr = acos( cz ) * degrad; - m_body.elevetr = 90.0 - m_body.zenetr; - refract(); - - // additional calculations for proper object sizing. - // orbit eccentricity - double e = 0.016709 - 1.151e-9 * daynumber; - // eccentric anomaly - double E = m_body.mnanom + e * degrad * sin(m_body.mnanom) * ( 1.0 + e * cos(m_body.mnanom) ); - double xv = cos(E) - e; - double yv = sqrt(1.0 - e*e) * sin(E); - m_body.distance = sqrt( xv*xv + yv*yv ); + m_body.zenetr = std::acos( cz ) * radtodeg; + m_body.elevetr = 90.0 - m_body.zenetr; + refract(); } void cSun::refract() { diff --git a/sun.h b/sun.h index e15abbbf..0834df73 100644 --- a/sun.h +++ b/sun.h @@ -66,23 +66,25 @@ protected: struct celestialbody { // main planet parameters - double dayang; // day angle (daynum*360/year-length) degrees - double mnlong; // mean longitude, degrees - double mnanom; // mean anomaly, degrees - double eclong; // ecliptic longitude, degrees. - double ecobli; // obliquity of ecliptic. - double declin; // declination--zenith angle of solar noon at equator, degrees NORTH. - double rascen; // right ascension, degrees - double hrang; // hour angle--hour of sun from solar noon, degrees WEST - double zenetr; // solar zenith angle, no atmospheric correction (= ETR) - double zenref; // solar zenith angle, deg. from zenith, refracted - double coszen; // cosine of refraction corrected solar zenith angle - double elevetr; // solar elevation, no atmospheric correction (= ETR) - double elevref; // solar elevation angle, deg. from horizon, refracted. - double distance; // distance from earth in AUs - double erv; // earth radius vector (multiplied to solar constant) - double etr; // extraterrestrial (top-of-atmosphere) W/sq m global horizontal solar irradiance - double etrn; // extraterrestrial (top-of-atmosphere) W/sq m direct normal solar irradiance + double dayang; // day angle (daynum*360/year-length) degrees + double phlong; // longitude of perihelion + double mnlong; // mean longitude, degrees + double mnanom; // mean anomaly, degrees + double tranom; // true anomaly, degrees + double eclong; // ecliptic longitude, degrees. + double oblecl; // obliquity of ecliptic. + double declin; // declination--zenith angle of solar noon at equator, degrees NORTH. + double rascen; // right ascension, degrees + double hrang; // hour angle--hour of sun from solar noon, degrees WEST + double zenetr; // solar zenith angle, no atmospheric correction (= ETR) + double zenref; // solar zenith angle, deg. from zenith, refracted + double coszen; // cosine of refraction corrected solar zenith angle + double elevetr; // solar elevation, no atmospheric correction (= ETR) + double elevref; // solar elevation angle, deg. from horizon, refracted. + double distance; // distance from earth in AUs + double erv; // earth radius vector (multiplied to solar constant) + double etr; // extraterrestrial (top-of-atmosphere) W/sq m global horizontal solar irradiance + double etrn; // extraterrestrial (top-of-atmosphere) W/sq m direct normal solar irradiance }; struct observer { // weather, time and position data in observer's location