List 1
/*
* Heat flow simulation model of multiple terrestrial geological columns
* Single 500 m ultrafine snow layer at 47.5 degrees N
*
* Source file: Glaciation201.java
*
* author: C F Davis
* version: 05 Oct 2019
*/
package glaciation;
import java.awt.*;
import java.util.*;
import java.awt.event.*;
import java.lang.Math;
import java.text.SimpleDateFormat;
import javax.swing.JTable;
import java.io.BufferedReader;
import java.io.FileReader;
import java.io.IOException;
import java.lang.Object;
public class Glaciation extends java.applet.Applet implements Runnable {
private Thread testThread = null;
EventHandler eventHandler;
int count, xmax, ymax;
Image offImage;
Graphics offGraphics;
Color colorpalette[] = new Color[15];
int xpoints[] = new int[4];
int ypoints[] = new int[4];
int x1, x2, y1, y2, paletteindex;
int p1, p2, p3, p4;
boolean go = true;
ItemGraphCanvas canvas1;
Image fieldImage;
Graphics fieldGraphics;
LayerDisplay layerDisplay;
Image layerImage;
Graphics layerGraphics;
Image textImage;
Graphics textGraphics;
int layerMargin;
boolean printManager = false;
int mode = 1;
int runmode = 1;
boolean clearScreen = true;
double timeOfDay = 0.0; // time of day 0 - 23
double hoursPerDay = 10.0;
double hourstimestep = 0.5; // time step interval (hrs)
int day = 0;
int outcount = 0;
// double conductionHeatGain = 0.0;
// double subtemp[][] = new double[2][500]; // subsurface temperatures
double tcount, tsum, tmean;
int selector = 0;
int opcontrol = 2;
// asteroid additions
// GeologicalColumn asteroid;
PseudoMap map;
GeologicalColumn column[] = new GeologicalColumn[20];
double secondsPerYear = 31540000.0;
double scale;
double screenWidth = 400;
Dimension appletSize;
double year, century, millennium;
long epochCZYN = 66043000; // nominal current epoch
long currentCZYN = 66043000; // Cenozoic Tear Bumber 0 is 66 million years before 1 Jan 2000 (Chicxulub crater)
long eemianCZYN = 66043000 - 100000; // Emian interglacial
long longCZYN = 66043000;
long startCZYN = eemianCZYN;
double initialRadius = 6371000;
double maxRadius = 6371000;
double initialLayerHeight = 17;
double fractionGranite = 1.0;
double heatGranite = 1E-9;
long tick = 0;
long longPeriodIndex = 0;
double previousGlaciationTotal = 0;
double glaciationTotal = 0; // counter
double glaciationPeriod = 0; // counter
double glaciationMaxPeriod = 1000; // period over which to measure glaciation
double elapsedSecs;
double glaciationPercent = 0; // % glaciation over glaciationMaxPeriod
boolean rollingGraphUpdated = false;
int rollingGraphWidth;
int rollingGraphHeight;
int rollingGraphStreams;
double longDay = 1.0 / 365.0; // there are 365 longDays in a year, 365.25 days in a year
int timeStepLongDays = 1; // timestep is an integer number of long days : valid numbers 1. 5, 73
// double timeStepYears = 1.0 / 365.25;
double timeStepYears = timeStepLongDays * (double)longDay;
// String csvFilename = "IcedHotAsteroid.csv";
// String csvFilename = "Earth1al5.csv";
/// String csvFilename = "Earth_20pcCC_45c.hcm";
// String csvFilename = "Ground1.csv";
// String csvFilename = "Earth_20cc_solar2_80.csv";
// String csvFilename = "E9.hcm";
// String csvFilename = "Earth_double_air_mass.hcm";
String csvFilename = "y120k.hcm";
// String csvFilename = "y65million.hcm";
double stime;
int degreesPerLatitude = 5;
int nLatitudes = 90 / degreesPerLatitude;
double viewLatitude = 47.5;
int nviewLatitude;
int solar = 2;
Milankovitch mil;
boolean milankovitch = false;
// double cloudCover = 0.20; // http://www.climate4you.com/ClimateAndClouds.htm
double cloudCover = 0.20;
boolean waterRunoff = true;
int icedropCounter = 0; // countdown starts when ice is dropped on surface
int icedropLayer = 67; // layer temperature to record after ice drop
double pauseDate = 449000;
boolean constantAddLayerDepther = false; // maintain added ice/snow layer depth
double addLayerMaterial = 3.05; // added layer material
// double addLayerMaterial = 11.0; // added layer material
double addLayerDepth = 20.0; // added layer depth
boolean relativeTemperatures = false; // relative temperature = current temperature - initial temperature
boolean heatEntireLayers = true; // true = changee layer pseudotemperature during phase change, false to melt off water
boolean varyIceCharacter = false; // true to vary ice conductivity with temperature
boolean varyAirEmissivity = false;
double addedIceMass = 0;
double addedIceHeatContentPerKg = 0;
double H20LayerHeatContent;
double H20LayerMassxSpht;;
double H20LayerMassxLatht;
double H20LayerHeatAdded;
double tempvar1 = 88.0; // any value other than 3.0 to activate old phase handling code
double snowHeatAdded = 0;
boolean albedoTest = false; // albedo test stepss thin surface layer albedo from 0 to 1
double meanGlobalTemperatureK;
double usefulVariable1 = 0; // variable that can be used in any way desired
double usefulVariable2 = 0; // variable that can be used in any way desired
DailyValues ocean[] = new DailyValues[ 2 * nLatitudes ]; // pre-calculated ocean daily air temperatures 36 latitudes
DailyValues land[] = new DailyValues[ 2 * nLatitudes ]; // pre-calculated ocean daily air temperatures 36 latitudes
DailyValues snow[] = new DailyValues[ 2 * nLatitudes ]; // pre-calculated ocean daily air temperatures 36 latitudes
DailyValues combined[] = new DailyValues[ 2 * nLatitudes ]; // combined daily air temperatures 36 latitudes
DailyValues active[] = new DailyValues[ 2 * nLatitudes ]; // active daily air temperatures used in surface temperature calculation
DailyValues global = new DailyValues(); // pre-calculated global daily mean air temperature
int nStartLatitude = 0;
int nStopLatitude = 17;
double annualGlobalMeanTemperature;
long suspendCalculationDate = 7000000; // debug aid: when not in use use value of 1000 million
boolean snowing = false;
long delayedStart = 12000;
int rollingMaxValue, rollingMinValue, rollingCalendarStep, rollingCalendarDisplayInterval;
double topLayerDepth = 40000;
int nband[] = new int[10];
boolean removeFlag = false;
boolean flag1 = false;
int videoWindow[][] = {
{ 360, 240, },
{ 480, 360, },
{ 640, 480, },
{ 1280, 720, },
{ 1920, 1080, }
};
public void init() {
stime = System.currentTimeMillis();
this.enableEvents( AWTEvent.MOUSE_EVENT_MASK );
double var = 0.4 / ( 10E3 * this.secondsPerYear);
System.out.println( "var " + var );
MaterialPhase.phase( 500, 100e3 );
// setSize( (int)screenWidth, (int)screenWidth );
int nScreen = 2;
screenWidth = videoWindow[nScreen][0];
setSize( videoWindow[nScreen][0], videoWindow[nScreen][1] );
appletSize = this.getSize();
xmax = appletSize.width; // graphics width
ymax = appletSize.height; // graphics height
setBackground(Color.white);
count=0;
year = 0;
century = 0;
millennium = 0;
nviewLatitude = (int)( ( this.viewLatitude - 2.5 ) / 5.0 );
// nStartLatitude = nviewLatitude;
// nStopLatitude = nviewLatitude;
System.out.println( "nviewLatitude " + nviewLatitude );
scale = videoWindow[nScreen][0] / ( maxRadius * 2.5 );
offImage = createImage( videoWindow[nScreen][0], videoWindow[nScreen][0] );
offGraphics = offImage.getGraphics();
layerDisplay = new LayerDisplay( this );
// set up textImage for text display
textImage = createImage( ( 3 * videoWindow[nScreen][0] / 4 ), ( videoWindow[nScreen][0] / 2 ) );
textGraphics = textImage.getGraphics();
map = new PseudoMap();
for ( int n=0; n<nLatitudes*2; n++ ) {
ocean[n] = new DailyValues();
land[n] = new DailyValues();
snow[n] = new DailyValues();
combined[n] = new DailyValues();
active[n] = new DailyValues();
}
// glaciate pseudoMap above viewLatitude
map.glaciateLandPseudoMap( this.nviewLatitude );
map.printPseudoMap();
// layer display setup
layerMargin = videoWindow[nScreen][0] / 4;
layerImage = createImage( layerMargin, videoWindow[nScreen][0] );
layerGraphics = layerImage.getGraphics();
layerDisplay.setDisplaySize( layerMargin, videoWindow[nScreen][1] );
System.out.println( "layerMargin " + layerMargin );
// rolling graph display setup
rollingGraphWidth = videoWindow[nScreen][0] * 3 / 4;
rollingGraphHeight = videoWindow[nScreen][1] / 2;
canvas1 = new ItemGraphCanvas( this, offGraphics, rollingGraphWidth, rollingGraphHeight, 5, true );
fieldImage = createImage( rollingGraphWidth, rollingGraphHeight );
fieldGraphics = fieldImage.getGraphics();
fieldGraphics.setColor(Color.white);
fieldGraphics.fillRect(0, 0, rollingGraphWidth, rollingGraphHeight );
canvas1.calendarLabel = "e5";
canvas1.calendarStep = 1000; // years per pixel
/* get displaymode */
try {
mode = Integer.parseInt( getParameter("MODE" ) );
} catch( Exception e ) {
mode = 1;
}
colorpalette[0] = new Color(000, 000, 000); // black
colorpalette[1] = new Color( 64, 64, 64); // light grey
colorpalette[2] = new Color(128, 128, 128); // mid gray
colorpalette[3] = new Color(192, 192, 192); // darke gray
colorpalette[4] = new Color(255, 255, 255); // white
colorpalette[5] = Color.RED;
colorpalette[6] = Color.ORANGE;
colorpalette[7] = Color.YELLOW;
colorpalette[8] = Color.CYAN;
colorpalette[9] = new Color( 85, 52, 52); // brown?
colorpalette[10] = new Color(135,206, 250); // sky blue
colorpalette[11] = new Color(240,240, 240); // steam
colorpalette[12] = new Color( 0, 0, 255); // blue
colorpalette[13] = new Color( 0, 0, 255); // blue
colorpalette[14] = new Color( 0, 0, 255); // blue
// asteroid = new GeologicalColumn( this, csvFilename, viewLatitude, 0 );
// asteroid = new GeologicalColumn( this, 0 );
// asteroid.setRadioactiveHeatGeneration( 2, mantle( 6371000, 3400000, 4650, 9990, 20E12 ) );
for ( int c=0; c<=17; c++ ) {
column[c] = new GeologicalColumn( this, csvFilename, (c*5)+2.5, 1 );
column[c].columnActiveCalculation = true;
}
// System.out.println( "Terrain albedo " + column[ nviewLatitude ].terrainAlbedo( 500, 1000, 0.4, 200, 0.8 ) );
if ( timeStepLongDays > 1 ) preliminaries();
// initialise graphics
this.pauseDate = 450000;
this.rollingMaxValue = 500;
this.rollingMinValue = 0;
this.rollingCalendarStep = 1000;
canvas1.calendarStep = this.rollingCalendarStep;
this.rollingCalendarDisplayInterval = 100000;
canvas1.calendarDisplayInterval = this.rollingCalendarDisplayInterval;
canvas1.calendarLabel = "e5";
this.topLayerDepth = 200000;
layerDisplay.topLayerDepth = this.topLayerDepth;
// initialise geological columns for runtime
this.timeStepLongDays = 1;
setTimeStepYears( this.timeStepLongDays );
this.viewLatitude = 47.5;
this.nviewLatitude = (int)( ( this.viewLatitude - 2.5 ) / 5.0 );
nStartLatitude = this.nviewLatitude;
nStopLatitude = this.nviewLatitude;
for ( int c=0; c<=17; c++ ) {
column[c] = new GeologicalColumn( this, csvFilename, (c*5)+2.5, 1 );
if ( c >= nStartLatitude && c <= nStopLatitude ) {
column[c].columnActiveCalculation = true;
} else {
column[c].columnActiveCalculation = false;
}
// explicitly initialise experiment variables
column[c].airMixingFactor = 0;
column[c].snowStartCZYN = this.startCZYN + 12000; // start snow after 12000 years
column[c].snowStopCZYN = this.startCZYN + 12240;
column[c].snowType = 3.05;
column[c].snowLayerThickness = 500;
column[c].snowDepositionInterval = 350;
column[c].terrainHeight = 0;
column[c].terrainWidth = 600;
// surface rock bands for rollingGraph
nband[0] = column[c].surfaceRockBand;
nband[1] = nband[0] - 1;
nband[2] = nband[1] - 1;
}
System.out.println( "nviewLatitude " + nviewLatitude );
eventHandler = new EventHandler( this );
eventHandler.insertEvent( eemianCZYN );
eventHandler.insertEvent( eemianCZYN + 1007 );
eventHandler.insertEvent( eemianCZYN + 15018 );
currentCZYN = eventHandler.eventQueue[0];
eventHandler.printEventQueue();
System.out.println( (long)(eemianCZYN + 15018) );
// milankovitchTest();
// mil = new Milankovitch( "INSOLN.LA2004.BTL.100.ASC" );
column[ nviewLatitude ].setRadioactiveHeatGeneration( 2, 2.1773010891400347E-8 );
// double r1 = 100.0;
// double r2 = 110.0;
// double v = column[ nviewLatitude ].sphereVolume(r2) - column[ nviewLatitude ].sphereVolume(r1);
// System.out.println( "r2 = " + column[ nviewLatitude ].layerRadius( v, column[ nviewLatitude ].sphereVolume(r1) ) );
tsum = 0.0; tcount = 0.0;
year = 0;
mil = new Milankovitch( this, "milankovitch250kyrs.csv", this.epochCZYN, 1000 );
}
void setTimeStepYears( int longDays ) {
if ( longDays==1 || longDays==5 || longDays==73 ) {
timeStepLongDays = longDays;
timeStepYears = timeStepLongDays * (double)longDay;
System.out.println( "timeStepYears " + timeStepYears );
}
}
public void start() {
if (testThread == null) {
testThread = new Thread(this, "Test1");
testThread.start();
}
}
public void processMouseEvent( MouseEvent e) {
if ( e.getID() == MouseEvent.MOUSE_ENTERED ) { go = false; }
else if ( e.getID() == MouseEvent.MOUSE_EXITED ) { go = true; }
else if ( e.getID() == MouseEvent.MOUSE_RELEASED ) {
offGraphics.setPaintMode();
offGraphics.setColor( Color.white );
offGraphics.fillRect( 0, 0, xmax, ymax );
x2 = e.getX();
y2 = e.getY();
p1 = 1 + x2 % 5;
p2 = 1 + y2 % 5;
p3 = 1 + Math.abs( x2-y2 ) % 5;
p4 = 1 + (x2 + y2) % 6;
go = true;
}
else super.processMouseEvent(e);
// System.out.println("The line number is " + new Exception().getStackTrace()[0].getLineNumber());
Object o = this;
Class c = o.getClass();
// System.out.println("class name is: " + c.getName());
// System.out.println("method name is: " + new Exception().getStackTrace()[0].getMethodName());
// System.out.println("calling method name is: " + Thread.currentThread().getStackTrace()[2].getMethodName());
// System.out.println("The class name is " + java.lang.Class.getSimpleName() );
}
public void run() {
runmode = mode;
Thread.currentThread().setPriority(Thread.MIN_PRIORITY);
Thread myThread = Thread.currentThread();
// preliminaries();
while (testThread == myThread) {
if ( go ) {
// if ( this.year > 15000 && this.year < 15001 ) this.setTimeStepYears( 5 );
mathEngine();
if ( opcontrol == 2 ) {
// clearScreen = true;
// repaint();
}
// if ( (currentCZYN - eemianCZYN) > 6.5942E7 ) System.out.println( year + " " + eventHandler.printAsteroidCZYN + " " + eventHandler.rollingGraphCZYN + " " + eventHandler.pauseEventCZYN );
runmode++; // change display graph mode
if ( runmode > 5) { runmode= 0; }
count++;
}
try {
Thread.sleep(0);
} catch (InterruptedException e){ }
}
}
// set all initial temperatures to current temperatures
void setInitialTemperatures() {
for ( int nlat=0; nlat<this.nLatitudes; nlat++ ) {
column[ nlat ].setInitialTemperatures();
}
}
void printAirTemperatureTables( int nlat ) {
int n = 0;
System.out.println( " day ocean land snow nlat=" + nlat );
for ( n=1; n<=365; n++ ) {
System.out.printf("% 5d ", n );
System.out.printf("%-10.2f", ocean[ nlat ].getDailyValue( n ) );
System.out.printf("%-10.2f", land[ nlat ].getDailyValue( n ) );
System.out.printf("%-10.2f", snow[ nlat ].getDailyValue( n ) );
System.out.println();
}
}
void mathEngine() {
boolean columnAlbedoChange = false;
double ytimestep = this.timeStepYears;
// do eventchecking first
eventHandler.checkEventQueue();
// 365 longDay year loop
for ( int dayIndex=this.timeStepLongDays; dayIndex<=365; dayIndex+= this.timeStepLongDays ) {
year += ytimestep;
century = year / 100.0;
millennium = year / 1000.0;
// if ( year > 100 && year < 101 ) System.out.println( year + " -- " + dayIndex );
// find new temperatures
for ( int nlat= 0; nlat<=17; nlat++ ) {
if ( column[ nlat ].columnActiveCalculation ) {
column[ nlat ].getNewLayerTemperatures( ( this.secondsPerYear * ytimestep ), this.timeStepLongDays, this.nviewLatitude, dayIndex );
if ( column[ nlat ].albedoChange ) {
columnAlbedoChange = true;
}
}
// at end of year, get new daily mean air temperature at latitude nlat
if ( dayIndex == 365 ) active[ nlat ].getMeanDailyValue();
}
tick += timeStepLongDays;
}
// if a column's albedo has changed, find new global air temperatures
if ( columnAlbedoChange ) findGlobalMeanAirTemperature();
/* code to remove top rock layers
if ( year > 12000 && !removeFlag ) {
// remove top two surface rock bands
column[ nviewLatitude ].removeLayer( column[ nviewLatitude ].surfaceRockBand );
column[ nviewLatitude ].reSynchronise();
column[ nviewLatitude ].removeLayer( column[ nviewLatitude ].surfaceRockBand );
column[ nviewLatitude ].reSynchronise();
this.nband[0] = column[ nviewLatitude ].surfaceRockBand;
removeFlag = true;
}
*/
currentCZYN++;
longCZYN++;
}
void rollingGraphOutput( double interval) {
String s;
int maxn = 5;
canvas1.day = (int)(interval);
// canvas1.calendarStep = (int)(interval*10); // years per pixel
// canvas1.calendarStep = 20; // years per pixel
double maxtemp = 500; //absolute values K
double mintemp = 0;
if ( relativeTemperatures ) {
maxtemp = 30;
mintemp = -20;
}
int n = this.nband[0];
// int n = column[ this.nviewLatitude ].surfaceRockBand;
double var = column[ this.nviewLatitude ].layerTable[ n ].averageTemperatureK.currentAverage;
if ( relativeTemperatures ) var = var - column[ this.nviewLatitude ].layerTable[ n ].initialTemperature;
canvas1.nexty[0] = var;
s = "" + 0;
canvas1.setStreamLabel(0, s);
canvas1.setStreamMaxMin( 0, maxtemp, mintemp );
// n--;
// var = column[ this.nviewLatitude ].layerTable[ n ].averageTemperatureK.currentAverage;
// if ( relativeTemperatures ) var = var - column[ this.nviewLatitude ].layerTable[ n ].initialTemperature;
n = this.nband[0] + 1; // air temperature
var = column[ this.nviewLatitude ].layerTable[ n ].currentTemperatureK;
// only plot snow temperature if snow/ice/water present
if ( var > 0 && column[ this.nviewLatitude ].layerTable[ n ].genericMaterialType() == 3 ) {
if ( flag1 == true ) {
canvas1.plot[1] = true;
} else {
canvas1.plot[1] = false; // don't plot surface rock temperature
flag1 = true;
}
} else {
canvas1.plot[1] = false;
flag1 = false;
}
canvas1.nexty[1] = var;
s = "" + 1;
canvas1.setStreamLabel(1, s);
canvas1.setStreamMaxMin( 1, maxtemp, mintemp );
// n--;
n = this.nband[0] - 1;
var = column[ this.nviewLatitude ].layerTable[ n ].averageTemperatureK.currentAverage;
if ( relativeTemperatures ) var = var - column[ this.nviewLatitude ].layerTable[ n ].initialTemperature;
canvas1.nexty[2] = var;
s = "" + 2;
canvas1.setStreamLabel(2, s);
canvas1.setStreamMaxMin( 2, maxtemp, mintemp );
// n = column[ this.nviewLatitude ].airLowestBand;
// var = column[ this.nviewLatitude ].layerTable[ n ].averageTemperatureK.currentAverage;
var = active[ this.nviewLatitude ].meanDailyValue;
// var = column[ this.nviewLatitude ].airLowestBandTemperatureK;
if ( relativeTemperatures ) var = var - column[ this.nviewLatitude ].layerTable[ n ].initialTemperature;
canvas1.nexty[3] = var;
s = "" + 3;
canvas1.setStreamLabel(3, s);
canvas1.setStreamMaxMin( 3, maxtemp, mintemp );
// n = column[ this.nviewLatitude ].thinSurfaceBand;
// var = column[ this.nviewLatitude ].layerTable[ n ].averageTemperatureK.currentAverage;
var = column[ this.nviewLatitude ].iceWaterDepth;
// if ( relativeTemperatures ) var = var - column[ this.nviewLatitude ].layerTable[ n ].initialTemperature;
canvas1.nexty[4] = var;
s = "" + 4;
canvas1.setStreamLabel(4, s);
canvas1.setStreamMaxMin( 4, maxtemp, mintemp );
}
public void paint(Graphics g) {
update(g);
}
public void update( Graphics g ) {
int n, bradius, nlayers;
long lyear;
int xcentre, ycentre, xtopleft, ytopleft, hr;
String str;
double topRadiusMaterial[] = new double[20];
if ( rollingGraphUpdated ) {
canvas1.paint( fieldGraphics );
this.rollingGraphUpdated = false;
} else {
// System.out.println( "roling graph not updated ");
}
// extract date from graphicLayer array of layer data
n = 0;
int material = -1;
while ( column[ nviewLatitude ].graphicLayer[n][0] >= 0 ) {
if ( column[ nviewLatitude ].graphicLayer[n][1] != material ) {
material = (int)column[ nviewLatitude ].graphicLayer[n][1];
topRadiusMaterial[ material ] = column[ nviewLatitude ].graphicLayer[n][2];
}
n++;
}
nlayers = n;
offGraphics.setPaintMode();
if ( clearScreen ) {
offGraphics.setColor( Color.white );
offGraphics.fillRect( 0, 0, xmax, ymax/2 ); // only clear top half of offImage
}
lyear = (long)year;
// display screen text
textGraphics.setPaintMode();
textGraphics.setColor( Color.white );
textGraphics.fillRect( 0, 0, xmax, ymax );
textGraphics.setColor( Color.black );
int px = 15;
int hpx = 200;
elapsedSecs = (System.currentTimeMillis() - stime)/1000.0;
if ( !go ) { str = " PAUSED" ; } else { str = ""; }
textGraphics.drawString( "CZYN " + (float)(currentCZYN - eemianCZYN) + str, 2, px );
textGraphics.drawString( "Elapsed secs " + elapsedSecs , hpx, px );
px += 15;
textGraphics.drawString( "air temperature " + (float)active[ this.nviewLatitude ].meanDailyValue, 2, px );
textGraphics.drawString( "CZY/sec " + (float)( (currentCZYN - eemianCZYN)/elapsedSecs ), hpx, px );
px += 15;
// textGraphics.drawString( "Rock " + (int)column[ nviewLatitude ].surfaceRockBand + " hfr=" + (float)column[ nviewLatitude ].surfaceRockHeatFlowRate + " (" + (float)column[ nviewLatitude ].maxSurfaceRockHeatFlowRate + ")", 2, px );
textGraphics.drawString( "Rock " + (int)column[ nviewLatitude ].surfaceRockBand + " hfr=" + (float)column[ nviewLatitude ].annualMeanHeatflow, 2, px );
textGraphics.drawString( "DT " + (float)(timeStepYears * 365.25) + " days", hpx, px );
px += 15;
textGraphics.drawString( "IceWater Depth " + (float)column[ nviewLatitude ].iceWaterDepth, 2, px );
textGraphics.drawString( "Filename " + csvFilename, hpx, px );
px += 15;
textGraphics.drawString( "Latitude " + viewLatitude, hpx, px );
double var;
boolean showTemperatures = true;
// show 3 surface rock layers
for ( n=column[ nviewLatitude ].graphicSurfaceRockBand-5; n<column[ nviewLatitude ].graphicLayerCount; n++ ) {
if ( column[ nviewLatitude ].graphicLayer[n][0] >= 0 ) {
if ( showTemperatures ) {
var = column[ nviewLatitude ].graphicLayer[column[ nviewLatitude ].graphicLayerCount-n-1][4];
if ( this.relativeTemperatures ) var = column[ nviewLatitude ].graphicLayer[column[ nviewLatitude ].graphicLayerCount-n-1][4] - column[ nviewLatitude ].graphicLayer[column[ nviewLatitude ].graphicLayerCount-n-1][9];
} else {
// show conductive heat flows to layer above
var = column[ nviewLatitude ].graphicLayer[column[ nviewLatitude ].graphicLayerCount-n-1][6];
}
textGraphics.drawString( "" + (int)column[ nviewLatitude ].graphicLayer[column[ nviewLatitude ].graphicLayerCount-n-1][1] + " UNI" + (int)column[ nviewLatitude ].graphicLayer[column[ nviewLatitude ].graphicLayerCount-n-1][7] + " T" + (int)column[ nviewLatitude ].graphicLayer[column[ nviewLatitude ].graphicLayerCount-n-1][0] + "= " + (float)var, 2, px );
px += 15;
}
}
// paint the asteroid layers into layerImage
layerDisplay.paint( layerGraphics, topRadiusMaterial[ 2 ] );
// layerDisplay.paint( layerGraphics, column[ nviewLatitude ], topRadiusMaterial[ 2 ] );
// draw the layerimage into top left corner of offImage
offGraphics.drawImage( layerImage, 0, 0, this );
// draw the textImage offset into offImage
offGraphics.drawImage( textImage, layerMargin, 0, this );
// draw the rolling graph fieldImage offset into offImage
offGraphics.drawImage( fieldImage, layerMargin, ymax/2, this );
g.drawImage( offImage, 0, 0, this );
}
Color setColour( int index ) {
Color colour = Color.white;
if ( index == 1 ) colour = Color.black;
if ( index == 2 ) colour = Color.darkGray;
if ( index == 3 ) colour = Color.lightGray;
if ( index == 4 ) colour = Color.cyan;
return( colour );
}
/*
Color setColour( int index , double temperature ) {
Color colour = Color.white;
Material mt = new Material( index );
int cindex = 0;
if ( temperature > mt.phaseTemperature1 ) cindex = 1;
if ( temperature > mt.phaseTemperature2 ) cindex = 2;
colour = colorpalette[ mt.colourIndex[ cindex ] ];
// if ( index == 1 ) colour = Color.black;
// if ( index == 2 ) colour = Color.darkGray;
return( colour );
}
*/
Color setColour( int index, double phase ) {
Color colour = Color.white;
Material mt = new Material( column[ nviewLatitude ], index );
int cindex = 0;
if ( phase == 1 ) cindex = 1;
if ( phase == 2 ) cindex = 2;
colour = colorpalette[ mt.colourIndex[ cindex ] ];
return( colour );
}
int scaleRadius( double r, double jscale ) {
int pixelHeight = (int)( r * jscale );
return( pixelHeight );
}
public void stop() {
testThread = null;
}
/*** start preliminaries *******************************************************/
// prelinary simulation to find air teemperatures at different latitudes and albedos
// using the simulation model
void preliminaries() {
int indexDay;
double initialTimestep = this.timeStepYears;
double y, latitudeAlbedo;
int y1 = 33;
int y2 = 66;
int y3 = 500;
// force timestep of one day
timeStepYears = 1.0 / 365.25;;
// enable calculation of atmosphere and thin surface layers
for ( int nlat=0; nlat<this.nLatitudes; nlat++ ) {
column[ nlat ].airMixingFactor = 0;
column[ nlat ].setActiveCalculationLayers ( column[ nlat ].topLayerOfAtmosphere , column[ nlat ].surfaceRockBand );
}
System.out.println( "setcCalculationLayers " + column[ 0 ].topLayerOfAtmosphere + " to " + column[ 0 ].surfaceRockBand );
int someTimePeriod = (int)( 1.0 / this.timeStepYears );
for ( int years=0; years<y3; years++ ) {
System.out.println( "preliminary year " + years );
latitudeAlbedo = 0.1; // ocean albedo
if ( years >= y1 ) latitudeAlbedo = 0.4 ; // granite/forest albedo
if ( years >= y2 ) latitudeAlbedo = 0.8 ; // snow albedo
indexDay = 1;
for ( int period=1; period<=someTimePeriod; period++ ) {
year += this.timeStepYears; // time step monthly
century = year / 100.0;
millennium = year / 1000.0;
for ( int nlat=0; nlat<this.nLatitudes; nlat++ ) {
column[ nlat ].setThinSurfaceBandAlbedo( latitudeAlbedo );
// find new temperatures, with surface rock temperatures fixed
column[ nlat ].getNewLayerTemperatures( ( this.secondsPerYear * this.timeStepYears ), 1, nlat, indexDay );
// duplicate N hemisphere temperatures offset 6 months in S hemispher
if ( latitudeAlbedo == 0.1 ) {
this.ocean[nlat].setDailyValue( column[ nlat ].airLowestBandTemperatureK, indexDay );
this.ocean[nlat + this.nLatitudes].setDailySouthernHemisphereValue( column[ nlat ].airLowestBandTemperatureK, indexDay );
} else if ( latitudeAlbedo == 0.4 ) {
this.land[nlat].setDailyValue( column[ nlat ].airLowestBandTemperatureK, indexDay);
this.land[nlat + this.nLatitudes].setDailySouthernHemisphereValue( column[ nlat ].airLowestBandTemperatureK, indexDay);
} else if ( latitudeAlbedo == 0.8 ) {
this.snow[ nlat ].setDailyValue( column[ nlat ].airLowestBandTemperatureK, indexDay);
// if ( nlat == this.nviewLatitude ) System.out.println( indexDay + ": " + (float)column[ nlat ].airLowestBandTemperatureK + ", " + (float)column[ nlat ].activeAirTemperature );
this.snow[nlat + this.nLatitudes].setDailySouthernHemisphereValue( column[ nlat ].airLowestBandTemperatureK, indexDay);
// if ( nlat + this.nLatitudes == 35 && years > 98 ) System.out.println( indexDay + " " + this.snow[nlat + this.nLatitudes].getDailyValue( indexDay ) );
}
}
// if ( year >= 98.9 ) System.out.println( (float)year + " " + indexDay + " " + currentGlobalMeanAirTemperature( indexDay ) );
tick++;
indexDay++;
}
}
// restore ap timestep
this.timeStepYears = initialTimestep;
this.year = 0;
System.out.println( "Annual global Mean Air Temperature " + findGlobalMeanAirTemperature() );
printAirTemperatureTables( nviewLatitude );
}
// find global mean air temperature as area-weighted sum of air temperatures above ocean, land, and snow
double currentGlobalMeanAirTemperature( int dayIndex ) {
double tm = 0;
double fractionSphereSurfaceArea;
for ( int nlat=0; nlat<this.nLatitudes; nlat++ ) {
fractionSphereSurfaceArea = map.pseudoMap[nlat+1][3];
// NORTHER hemisphere
tm += ocean[ nlat ].getDailyValue( dayIndex ) * map.pseudoMap[nlat+1][0] * fractionSphereSurfaceArea; // ocean
tm += land[ nlat ].getDailyValue( dayIndex ) * map.pseudoMap[nlat+1][1] * fractionSphereSurfaceArea; // land
tm += snow[ nlat ].getDailyValue( dayIndex ) * map.pseudoMap[nlat+1][2] * fractionSphereSurfaceArea; // snow
// SOUTHERN hemisphere
tm += ocean[ nlat + this.nLatitudes ].getDailyValue( dayIndex ) * map.pseudoMap[nlat+1 + this.nLatitudes][0] * fractionSphereSurfaceArea; // ocean
tm += land[ nlat + this.nLatitudes ].getDailyValue( dayIndex ) * map.pseudoMap[nlat+1 + this.nLatitudes][1] * fractionSphereSurfaceArea; // ocean
tm += snow[ nlat + this.nLatitudes ].getDailyValue( dayIndex ) * map.pseudoMap[nlat+1 + this.nLatitudes][2] * fractionSphereSurfaceArea; // ocean
}
return( tm );
}
// using Preliminary run
// get annual mean temperature from daily mean temperatures
double findGlobalMeanAirTemperature() {
double meanAirTemperature;
double totalGlobalMeanTemperature = 0;
for ( int dayIndex=1; dayIndex<=365; dayIndex++ ) {
meanAirTemperature = currentGlobalMeanAirTemperature( dayIndex );
totalGlobalMeanTemperature += meanAirTemperature;
this.global.setDailyValue( meanAirTemperature, dayIndex );
// System.out.println( "day " + dayIndex + " global mean air temperature " + meanAirTemperature );
}
this.global.getMeanDailyValue();
totalGlobalMeanTemperature = totalGlobalMeanTemperature / 365.0;
this.annualGlobalMeanTemperature = totalGlobalMeanTemperature;
System.out.println( year + " new annualGlobalMeanTemperature " + this.annualGlobalMeanTemperature);
return( totalGlobalMeanTemperature );
}
// update PseudoMap with new snow area, and find new combided air temperature
double combinedAirTemperature( int nlat, double snowCoverageFraction, int dayIndex ) {
double ta = 0;
double landPlusSnow = 1.0 - map.pseudoMap[nlat+1][0];
double snowArea = snowCoverageFraction * landPlusSnow;
double landArea = landPlusSnow - snowArea;
// update map
map.pseudoMap[nlat+1][2] = snowArea;
map.pseudoMap[nlat+1][1] = landArea;
// find combined ait temperature ta, and store in combined[] table
ta += ocean[ nlat ].getDailyValue( dayIndex ) * map.pseudoMap[nlat+1][0]; // ocean
ta += land[ nlat ].getDailyValue( dayIndex ) * map.pseudoMap[nlat+1][1] ; // land
ta += snow[ nlat ].getDailyValue( dayIndex ) * map.pseudoMap[nlat+1][2] ; // snow
this.combined[ nlat ].setDailyValue( ta, dayIndex );
if ( dayIndex == 365 ) this.combined[ nlat ].getMeanDailyValue();
// if ( year > 11999.5 && year < 12000.5 ) System.out.println( (float)snowArea + " " + (float)landArea + " : " + (float)ta );
return( ta );
}
/** end preliminaries *******************************************************/
}
class GeologicalColumn {
Glaciation ap;
boolean columnActiveCalculation;
double columnLatitude;
int nLatitude;
double areaAtLatitude; // area of circular strip at column latitude
// double fractionAreaAtLatitude; // fraction of asteroid surface area of strip at column latitude
// double airTemperatureAtLatitude;
// double fractionLandAreaAtLatitude = 0.6; // 60% land in N hemisphere 0% in S hemisphere
SolarRadiation sr;
double radius;
double nominalRadius; // radius used in graphic display
int nlayers;
LinkedList<ConductiveLayer> llayer;
double initialTemperatureK;
int nIterations = 0;
double referenceRadius = 0.0; // earth radius
int surfaceBand; // surface band for solar irradiation
int thinSurfaceBand; // thin radiative surface band number
// DailyValues glaciated = new DailyValues(); // pre-calculated glaciated daily air temperatures at colum latitude
// DailyValues unglaciated = new DailyValues(); // .. .. unglaciated .. .. .. .. .. ..
double stefanBoltzmannConstant = 5.6699E-8;
double thinSurfaceBandTemperature; // thin radiative surface band temperature
double thinSurfaceRadius = 0;
int airLowestBand; // bottom layer of air in atmospher
double airLowestBandTemperature; // temperature of bottom layer of air
double airLowestBandTemperatureK; // relative .. .. .. ..
double surfaceBandTemperature;
double lastSurfaceTemperature;
double surfaceRockHeatFlowRate;
double maxSurfaceRockHeatFlowRate = 0;
double surfaceRockTemperature = 0;
double atmosphereTemperature;
int surfaceRockBand;
double surfaceHeatFlowRate = 0;
double year;
double graphicLayer[][] = new double[2000][10]; // date for graphics display
int graphicLayerCount = 0; // data for grapjic display
int graphicSurfaceRockBand = 0; // data for graphic display
boolean resynchronise = true;
double annualSolarInsolation = 1362 * 3.6E6; // annual solar radiation kwh-> J UK
double solarGain; // solar gain Watts
double albedo = 0.3; // mean Earth albedo
double previousAlbedo = albedo; // albedo on previous iteration
boolean albedoChange = true; // set true at initialisation, and when albedo changes
double secondsPerYear = 3.154e+7;
double secondsPerDay = 86400;
double iceDepth = 0;
double iceWaterDepth = 0;
double previousIceWaterDepth = 0;
boolean iceMeltingAtBase = false;
boolean printFlag = false;
double surfaceArea; // surface area of asteroid
double mass; // mass of asteroid (earth 5.972 × 10^24 kg)
double atmosphereMass; // (earth 5.1480×10^18 kg)
double topOfAtmosphere;
int topLayerOfAtmosphere;
double surfaceGravity = 9.81; // surface gravitaional acceleration (should be calculated)
double asteroidHeatContent;
double meanTimestepSolarGain; // J / timestep
boolean varyIceAlbedo = true;
boolean varyIceAlbedoReset = true;
double phaseChangeDate[] = new double[200];
int uniqueNumber = 1;
int removeLayerFlag = -1; // number of layer to remove
boolean flagPrintColumn = false;
boolean flagSetAirEmissivity = false;
double yearSetAirEmissivity = 22000; // year that emissivity is to be set
double setAirEmissivity = 0.909; // high emissivity/absorptivity to emulate global warming
double startOfGlaciation;
double endOfGlaciation;
boolean setTSLtempEqualToSRtemp = false; // set thin surface layer temperature equal to top surface rock temperature (zero surface heat loss)
boolean fixedIceHeatflow = false; // flag fixed heat flow rate into ice/snow from surface rock below
double fixedIceHeatflowRate = 0.060; // fixed heat flow rate into ice/snow W
double forestHeight = 0.0;
double forestAlbedo = 0.09;
double thinSurfaceLayerAlbedo = 0.09;
int test;
double degreesToRadians = Math.PI / 180;
int timer = 0;
double numAnnualMeanValues = 0;
// double annualMeanAirTemperature;
double annualMeanHeatflow;
// double annualMeanAirTemperatureSum = 0;
double annualMeanHeatflowSum = 0;
ConductiveLayer layerTable[] = new ConductiveLayer[400];
double airMixingFactor = 0.001; // air mixing factot in range 0 ( no mixing) to 1 (complete mixing)
double activeAirTemperature;
double horizontalProportionOfSnow;
double snowLayerThickness = 25; // column snow deposition
double snowDepositionInterval = 275;
double snowMaximumDepth = 3000;
double snowType = 3.05;
double snowTemperature;
long snowStartCZYN = 0;
long snowStopCZYN = 66043000 * 2;
double precipitstionCZYN ;
double terrainHeight = 300;
double terrainWidth = 600;
double terrainUnderlyingAlbedo = 0.4;
double terrainSnowAlbedo = 0.8;
GeologicalColumn() { }
GeologicalColumn( Glaciation a, String filename, double lat, int test ) {
this.ap = a;
this.test = test;
System.out.println();
System.out.println( "column latitude " + lat );
this.columnLatitude = lat;
referenceRadius = 6371500; // Earth-specific
nominalRadius = referenceRadius;
topOfAtmosphere = referenceRadius;
System.out.println( "topOfAtmosphere0 " + topOfAtmosphere );
areaAtLatitude = AreaAtLatitude( this.columnLatitude, 2.5 );
System.out.println( "area At " + (float)this.columnLatitude + " Latitude " + areaAtLatitude );
this.surfaceArea = this.sphereSurfaceArea( referenceRadius );
// this.fractionAreaAtLatitude = this.areaAtLatitude / this.surfaceArea;
// System.out.println( "this.fractionAreaAtLatitude " + this.fractionAreaAtLatitude );
if ( filename.contains( ".csv" ) ) csvReader( filename );
else if ( filename.contains( ".hcm" ) ) hcmReader( filename );
else System.out.println( "no file" );
reSynchronise();
System.out.println( "topOfAtmosphere1 " + topOfAtmosphere );
// adjustLayerTemperatures( 270.0/281.0 );
// adjustLayerTemperatures( 0.5 );
// setLayerPhases();
// fixBandTemperature( true, 64 );
// if ( this.columnLatitude == this.ap.viewLatitude ) printColumn( 0 );
// if ( this.columnLatitude == this.ap.viewLatitude ) printColumnRC( 0 );
recalculateRadii( 0 );
// if ( this.columnLatitude == this.ap.viewLatitude ) printColumn( 0 );
// if ( this.airMixingFactor > 0 ) setActiveCalculationLayers( this.thinSurfaceBand, this.thinSurfaceBand - 20 );
// precipitstionCZYN = (long)ap.eemianCZYN + 12000;
this.precipitstionCZYN = (long)ap.startCZYN + ap.delayedStart;
// ap.eventHandler.insertEvent( ap.eventHandler.precipitstionCZYN );
// setLayerTicks( 65, 100 );
/*
System.out.println( "addLayer");
addLayer( 72, 270.0, 200.0, 3.05);
recalculateRadii( 0 );
printcsvColumn( 0 );
*/
double SRtimestep = 1.0 / 365.25;
sr = new SolarRadiation( this.ap, this, 1367.0, lat, SRtimestep, 365.25, 0.016702, 0.409093, 1.796257 );
System.out.println( "mean annual solar power: " + this.sr.meanSolarPower );
this.areaAtLatitude = areaAtLatitude();
// if ( this.columnLatitude == this.ap.viewLatitude ) setGeothermalGradient( 0.060 );
System.out.println( "topOfAtmosphere2 " + topOfAtmosphere );
}
// read .csv (comma separated variables) file containing body locations, velocities.
// (need to set working directory and applet.policy )
double csvReader( String file ) {
double m, x, y, z, xx, yy, zz, v, r, jd, startJulianDate;
double nband, mt, rLayer, wLayer, hLayer, layerTemperature;
double alayer[][] = new double[150][10];
ConductiveLayer layer;
llayer = new LinkedList<ConductiveLayer>();
int b = 0;
int rowcount, fieldcount, clr;
boolean bodyok, filefound = false;
String full_name, nm, dataRow;
BufferedReader CSVFile;
String rs[] = new String[400];
startJulianDate = 0.0;
rowcount = 0;
boolean atmosphereFound = false;
this.nlayers = 0;
try {
CSVFile = new BufferedReader(new FileReader(file));
rowcount = 0;
try {
dataRow = CSVFile.readLine(); // Read first line.
} catch (IOException e) {
dataRow = new String("");
}
while (dataRow != null) {
bodyok = false;
filefound = true;
StringTokenizer st = new StringTokenizer(dataRow, ",");
fieldcount = 0;
while (st.hasMoreTokens()) {
String s = st.nextToken();
rs[ fieldcount] = s;
fieldcount++;
}
if (rowcount > 0) {
nband = Double.valueOf(rs[0]).doubleValue();
if ( nband >= 0 ) {
// read layer data
alayer[nlayers][0] = Double.valueOf(rs[0]).doubleValue();
alayer[nlayers][1] = Double.valueOf(rs[1]).doubleValue();
alayer[nlayers][2] = Double.valueOf(rs[2]).doubleValue();
alayer[nlayers][3] = Double.valueOf(rs[3]).doubleValue();
alayer[nlayers][4] = Double.valueOf(rs[4]).doubleValue();
alayer[nlayers][5] = Double.valueOf(rs[5]).doubleValue();
alayer[nlayers][6] = Double.valueOf(rs[6]).doubleValue();
// single layer atmosphere is lowest air layer
if ( alayer[nlayers][1] == 5 ) {
topOfAtmosphere += alayer[nlayers][2];
topLayerOfAtmosphere = (int)nband;
// single layer atmosphere is first air mass found
if ( !atmosphereFound ) {
atmosphereFound = true;
}
// System.out.println( "air " + nband );
}
this.nlayers++;
} else if ( nband == -2 ) {
// read asteroid data
// this.referenceRadius = Double.valueOf(rs[1]).doubleValue();
}
}
rowcount++;
try {
dataRow = CSVFile.readLine(); // Read next line.
} catch (IOException e) {
dataRow = new String("");
}
}
// Close the file once all data has been read.
try {
CSVFile.close();
} catch (IOException e) {
}
} catch (IOException e) {
filefound = false;
System.out.println( file + " CSV file not found");
}
if ( filefound ) {
// construct asteroid column from centre outward to surface
int nl = nlayers;
for ( int i=nl-1; i>=0; i-- ) {
if ( alayer[i][0] >= 0 ) {
// alayer[i][6] = 5000;
// radius, width, height, temperature, material, tickDuration, W/m^3
layer = new ConductiveLayer( this, alayer[i][2], alayer[i][4], alayer[i][3], alayer[i][6], alayer[i][1], alayer[i][5], 0 );
layer.currentLayerPressure = layer.currentLayerPressure * 0.50;
llayer.addFirst( layer );
if ( (int)alayer[i][1] == 0 ) {
this.thinSurfaceRadius = alayer[i][2];
this.referenceRadius = this.thinSurfaceRadius;
System.out.println( "thinSurfaceRadius " + thinSurfaceRadius );
}
}
}
// End the printout with a blank line.
System.out.println("CSV File complete + " + rowcount + " rows " );
}
return( startJulianDate );
}
// read .hcm (comma separated variables) file containing layer heat contents and masses
// (need to set working directory and applet.policy )
double hcmReader( String file ) {
double nband, startJulianDate;
// double nband, mt, rLayer, wLayer, hLayer, layerTemperature;
double alayer[][] = new double[150][10];
ConductiveLayer layer;
llayer = new LinkedList<ConductiveLayer>();
int b = 0;
int rowcount, fieldcount, clr;
boolean bodyok, filefound = false;
String full_name, nm, dataRow;
BufferedReader CSVFile;
String rs[] = new String[400];
startJulianDate = 0.0;
rowcount = 0;
boolean atmosphereFound = false;
this.nlayers = 0;
System.out.println( file + " hcm file reader");
try {
CSVFile = new BufferedReader(new FileReader(file));
rowcount = 0;
try {
dataRow = CSVFile.readLine(); // Read first line.
} catch (IOException e) {
dataRow = new String("");
}
while (dataRow != null) {
bodyok = false;
filefound = true;
StringTokenizer st = new StringTokenizer(dataRow, ",");
fieldcount = 0;
while (st.hasMoreTokens()) {
String s = st.nextToken();
rs[ fieldcount] = s;
fieldcount++;
}
if (rowcount > 0) {
nband = Double.valueOf(rs[0]).doubleValue();
if ( nband >= 0 ) {
// read layer data
alayer[nlayers][0] = Double.valueOf(rs[0]).doubleValue(); // layer nBand
alayer[nlayers][1] = Double.valueOf(rs[1]).doubleValue(); // layer materialType
alayer[nlayers][2] = Double.valueOf(rs[2]).doubleValue(); // layer radius
alayer[nlayers][3] = Double.valueOf(rs[3]).doubleValue(); // layer height
alayer[nlayers][4] = Double.valueOf(rs[4]).doubleValue(); // layer mass
alayer[nlayers][5] = Double.valueOf(rs[5]).doubleValue(); // layer tick duration (secs)
alayer[nlayers][6] = Double.valueOf(rs[6]).doubleValue(); // layer heat content per Kg
alayer[nlayers][7] = Double.valueOf(rs[7]).doubleValue(); // layer temperature
// single layer atmosphere is lowest air layer
if ( alayer[nlayers][1] == 5 ) {
topOfAtmosphere += alayer[nlayers][2];
topLayerOfAtmosphere = (int)nband;
// single layer atmosphere is first air mass found
if ( !atmosphereFound ) {
atmosphereFound = true;
}
// System.out.println( "air " + nband );
}
this.nlayers++;
} else if ( nband == -2 ) {
// read asteroid data
// this.referenceRadius = Double.valueOf(rs[1]).doubleValue();
}
}
rowcount++;
try {
dataRow = CSVFile.readLine(); // Read next line.
} catch (IOException e) {
dataRow = new String("");
}
}
// Close the file once all data has been read.
try {
CSVFile.close();
} catch (IOException e) {
}
} catch (IOException e) {
filefound = false;
System.out.println( file + " hcm file not found");
}
if ( filefound ) {
// construct asteroid column from centre outward to surface
int nl = nlayers;
double r = 0;
for ( int i=nl-1; i>=0; i-- ) {
if ( alayer[i][0] >= 0 ) {
// nBand material radius height mass tickduration heat content temperature
layer = new ConductiveLayer( this, (int)alayer[i][0], alayer[i][1], r, alayer[i][3], alayer[i][4], alayer[i][5], alayer[i][6], alayer[i][7] );
layer.currentLayerPressure = layer.currentLayerPressure * 0.50;
r = layer.radius; // use current layer radius to calculate next layer radius
llayer.addFirst( layer );
if ( (int)alayer[i][1] == 0 ) {
this.thinSurfaceRadius = alayer[i][2];
this.referenceRadius = this.thinSurfaceRadius;
// System.out.println( "thinSurfaceRadius " + thinSurfaceRadius );
}
}
}
// End the printout with a blank line.
System.out.println("hcm File complete + " + rowcount + " rows " );
}
return( startJulianDate );
}
// find all condutctive heat flows
// using link list iterator
void getNewLayerTemperatures( double timeStepSeconds, int timeStepLongDays, int nlat, int dayIndex ) {
int removeLayer = -1;
int n, bmax, bmin, i;
boolean flag = false;
boolean continuation = true;
double heatflowrate, heatgain, dtemp, tmax, tmin, hfr, effectiveTemperature;
ConductiveLayer lyr[] = new ConductiveLayer[3];
int counter = 0;
ListIterator<ConductiveLayer> itr = null;
boolean timedOut = false;
this.nLatitude = nlat; // nLatitude should be set at initialisation
// if ( ap.year == 109.82067077339573 ) System.out.println( nlatitude + " lat in yr 109.82067077339573" );
// get the top 3 layers
// layer 0 is a dummy floating layer
heatflowrate = 0.0;
itr = llayer.listIterator();
if ( itr.hasNext() ) {
lyr[0] = itr.next();
counter++;
}
// layer 1 is a layer of air
if ( itr.hasNext() ) {
lyr[1] = itr.next();
counter++;
}
// layer 2 is either a layer of water or ice or rock
// lyr[] contains layer[n-1], layer[n], layer[n+1] for all layers
// While loop continues down from top of atmosphere until continuation is false,
// which happens when the temperature of the layer is fixed.
while(itr.hasNext() && continuation ){
lyr[2] = itr.next();
counter++;
// here come the heat flow calculations
heatflowrate = 0.0;
lyr[0].layerResistance();
lyr[1].layerResistance();
lyr[2].layerResistance();
lyr[1].updated = false;
// if ( lyr[1].nBand == this.thinSurfaceBand ) System.out.println( lyr[1].nBand + " thinSurfaceBand " + lyr[1].currentTemperatureK + " albedo " + lyr[1].albedo );
// if ( lyr[1].nBand == this.surfaceRockBand ) System.out.println( lyr[1].nBand + " thinSurfaceBand " + lyr[1].currentTemperatureK );
if ( lyr[1].genericMaterialType() == 0 && lyr[1].activeCalculation ) {
// thin surface layer
if ( ap.solar == 3 ) {
// Milankovitch daily solar power W/m^2
solarGain = ( sr.getTotalSolarGain( dayIndex-timeStepLongDays, dayIndex ) ) / timeStepSeconds; // assumes daily solar gain
} else if ( ap.solar == 2 ) {
solarGain = ( sr.getTotalSolarGain( dayIndex-timeStepLongDays, dayIndex ) ) / timeStepSeconds; // assumes daily solar gain
// if ( ap.year > 0 && ap.year < 2 ) System.out.println( ap.year + " dayIndex " + dayIndex + " timeStepLongDays " + timeStepLongDays + " solarGain " + solarGain );
} else if ( ap.solar == 1 ) {
solarGain = sr.meanSolarPower; // mean solar power (Watts) over year
} else {
solarGain = 0;
}
solarGain = (1.0 - ap.cloudCover) * solarGain; // solar gain is reduced by cloud cover
// if thin surface layer albedo < 0, use albedo of layer below
if ( lyr[1].albedo < 0 ) {
this.albedo = Math.abs( lyr[2].albedo );
// height width albedo snow depth snow albedo
this.albedo = this.terrainAlbedo( this.terrainHeight, this.terrainWidth, this.terrainUnderlyingAlbedo, this.iceWaterDepth, this.terrainSnowAlbedo, dayIndex );
} else {
this.albedo = Math.abs( lyr[1].albedo );
}
// test terrain albedo
// height width albedo snow depth snow albedo
// this.albedo = this.terrainAlbedo( 0, 600, 0.4, this.iceWaterDepth, 0.8, dayIndex );
this.thinSurfaceLayerAlbedo = this.albedo;
double ts = surfaceTemperature( lyr[0], lyr[1], lyr[2], dayIndex, timeStepLongDays ); // fourth root solution
if ( !lyr[1].fixTemperature ) {
lyr[1].layerHeatContentPerKg = ts * lyr[1].layerSpecificHeat;
lyr[1].layerHeatContent = lyr[1].layerHeatContentPerKg * lyr[1].mass;
this.lastSurfaceTemperature = ts;
lyr[1].nextTemperatureK = ts;
// System.out.println( "thin surface layer phase " + lyr[1].phase );
} else {
continuation = false;
}
// if ( ap.year >= 109.82067077339573 && ap.year < 109.99 ) System.out.println( ts + " ts " + this.columnLatitude + " columnLatitude " + ap.viewLatitude + " viewLatitude in yr " + ap.year );
} else if ( lyr[1].genericMaterialType() == 5 && lyr[1].activeCalculation ) {
// multi-layer atmosphere
// radiative heat exchanges from atmospheric air
heatflowrate = 0;
// heat radiated from surface above is all absorbed
heatflowrate += wattsRadiated( lyr[0].currentTemperatureK, lyr[0].layerEmissivity );
// heat radiated from surface below is not all absorbed. Some is radiated into space
heatflowrate += lyr[1].layerEmissivity * wattsRadiated( lyr[2].currentTemperatureK, lyr[2].layerEmissivity );
// heat radiated to surface below
heatflowrate -= wattsRadiated( lyr[1].currentTemperatureK, lyr[1].layerEmissivity );
// heat radiated to surface above
heatflowrate -= wattsRadiated( lyr[1].currentTemperatureK, lyr[1].layerEmissivity );
// find heat gain over timestep
// heatgain = heatflowrate * interval;
heatgain = heatflowrate * timeStepSeconds;
// find change in temperature due to heat gain.
// * 5.0 is a fix to make current 5 atmosphere layers less responsive and stable
dtemp = heatgain / ( lyr[1].mass * lyr[1].layerSpecificHeat * 1.0 );
if ( !lyr[1].fixTemperature ) {
lyr[1].nextTemperatureK = lyr[1].currentTemperatureK + dtemp;
} else {
continuation = false;
}
// System.out.println( "nband " + lyr[1].nBand + ": heatgain " + heatgain + " dtemp " + dtemp + " lyr[1].nextTemperatureK " + lyr[1].nextTemperatureK );
/*
timer++;
if ( timer >= 10000 ) {
System.out.println( ap.year + " this.nBand " + lyr[1].nBand );
timer = 0;
}
*/
} else if ( lyr[1].genericMaterialType() > 0 && lyr[1].activeCalculation ) {
// find if any snow above surface rocks is melthing
if ( lyr[1].nBand == this.surfaceRockBand ) {
/*
if ( lyr[0].genericMaterialType() == 3 && lyr[0].currentTemperatureK >= 273.0 ) {
this.iceMeltingAtBase = true;
} else {
this.iceMeltingAtBase = false;
}
*/
if ( lyr[1].currentTemperatureK > 273.0 ) {
this.iceMeltingAtBase = true;
} else {
this.iceMeltingAtBase = false;
}
}
// while layer phase change is in process, currentTemperatureK is fixed, and
// phaseTemperature is used instead until phase change is complete
effectiveTemperature = lyr[1].currentTemperatureK;
// conductive heat transfer
// heat flow rate rate from higher layer 0 to layer 1
if ( fixedIceHeatflow && ( lyr[0].genericMaterialType() == 3 ) && ( lyr[1].nBand == this.surfaceRockBand ) ) {
// here if layer is surface rock band and layer above is snow/ice
heatflowrate -= fixedIceHeatflowRate; // assumes upward heat flow
} else {
// here for normal heat flow
heatflowrate += conductiveHeatFlow( lyr[1], lyr[0].layerResistance, lyr[1].layerResistance, lyr[0].currentTemperatureK, effectiveTemperature );
}
// if ( lyr[1].nBand == 28 ) System.out.print( (float)lyr[0].currentTemperatureK + " 0->1 " + (float)heatflowrate );
lyr[1].upwardHeatFlowRate = heatflowrate;
if ( Math.abs( lyr[1].upwardHeatFlowRate ) > Math.abs( lyr[1].maxUpwardHeatFlowRate ) ) lyr[1].maxUpwardHeatFlowRate = lyr[1].upwardHeatFlowRate;
if ( lyr[1].nBand == this.surfaceRockBand ) {
this.surfaceRockHeatFlowRate = heatflowrate;
this.maxSurfaceRockHeatFlowRate = lyr[1].maxUpwardHeatFlowRate;
}
// conductive heat flow rate from lower layer layer 2 to layer 1
if ( fixedIceHeatflow && ( lyr[1].genericMaterialType() == 3 ) && ( lyr[2].nBand == this.surfaceRockBand ) ) {
// here if fixed heat flow rate into ice/snow from surface rock beneth
heatflowrate += fixedIceHeatflowRate; // assumes upward heat flow
} else {
// here for normal heat flow
heatflowrate += conductiveHeatFlow( lyr[1], lyr[2].layerResistance, lyr[1].layerResistance, lyr[2].currentTemperatureK, effectiveTemperature );
}
heatflowrate += lyr[1].heatProduced; // radioactive heat production in layer
// find heat gain over timestep
heatgain = heatflowrate * timeStepSeconds;
// heat a layer if ( lyr[1].nBand == this.surfaceRockBand-4 ) heatgain = Math.abs( heatgain );
lyr[1].layerHeatAdded += heatgain;
lyr[1].layerHeatContentPerKg += heatgain / lyr[1].mass;
lyr[1].layerHeatContent += heatgain;
// find change in temperature due to heat gain.
// set the next layer temperature that will replace current layer tmperature
lyr[1].nextTemperatureK = lyr[1].pt.temperature( lyr[1].layerHeatContentPerKg );
}
lyr[1].updated = true;
// roll layers up
lyr[0] = lyr[1];
lyr[1] = lyr[2];
}
// now copy next temperatures to current temperatures
// and perform other housekeeping tasks
int surfaceLayerIndex;
double asteroidRadius = 0;
asteroidHeatContent = 0;
counter = 0;
itr=llayer.listIterator();
while(itr.hasNext()){
lyr[0] = itr.next();
if ( lyr[0].updated ) {
// if ( lyr[0].nBand == 71 ) System.out.println( "l71 T " + lyr[0].currentTemperatureK );
// set current temperature of layer
lyr[0].previousTemperaureK = lyr[0].currentTemperatureK;
lyr[0].currentTemperatureK = lyr[0].nextTemperatureK;
lyr[0].averageTemperatureK.addNewValue( lyr[0].currentTemperatureK );
// if ( lyr[0].nBand == this.surfaceRockBand ) System.out.println( "this.surfaceRockBand " + lyr[0].averageTemperatureK.currentAverage );
/*
// store average surface rock temperatures ( topmost = 0)
surfaceLayerIndex = this.surfaceRockBand - lyr[0].nBand;
if ( surfaceLayerIndex >= 0 && surfaceLayerIndex < 10 ) {
this.surfaceTemperature[ surfaceLayerIndex ] = lyr[0].averageTemperatureK.currentAverage;
}
*/
// if phase has changed
if ( lyr[0].phase != (int)lyr[0].nextPhase ) {
lyr[0].phase = (int)lyr[0].nextPhase;
double mtrl = lyr[0].materialType;
mtrl = (int)mtrl + ( lyr[0].phase / 10.0 ); // 3.0 is ice, 3.1 is water, 3.2 is steam
lyr[0].materialType = mtrl;
lyr[0].material.setMaterialCharacteristics( mtrl, 1 );
// System.out.println( (float)ap.year + " bnnd " + lyr[0].nBand + " materialType becomes " + lyr[0].materialType );
// flag liquid or vapour water layer to be removed
if ( (int)mtrl == 3 ) {
if ( ap.waterRunoff && lyr[0].materialType >- 3.1 ) {
removeLayerFlag = lyr[0].nBand;
}
}
}
// find layer heat content, and total asteroid heat content
// asteroidHeatContent += lyr[0].heatContent();
// ice thermal characteristics vary with temperature (not during phase changes)
if ( ap.varyIceCharacter && lyr[0].material.isIce( lyr[0].currentTemperatureK) && !lyr[0].phaseChange ) {
i = lyr[0].material.indexIceTemperature( lyr[0].currentTemperatureK );
lyr[0].layerDensity = lyr[0].material.iceCharacteristics[i][1];
lyr[0].layerConductivity = lyr[0].material.iceCharacteristics[i][2];
lyr[0].layerSpecificHeat = lyr[0].material.iceCharacteristics[i][3] * 1000.0;
lyr[0].layerResistance = lyr[0].layerResistance();
}
// find temperature of lowest atmosphere layer (last air temperature found)
if ( lyr[0].material.isAir( lyr[0].currentTemperatureK ) ) {
this.atmosphereTemperature = lyr[0].currentTemperatureK;
// if ( ap.relativeTemperatures ) this.atmosphereTemperature = lyr[0].currentTemperatureK - lyr[0].initialTemperature ;
}
// set surface rock temperature.
if ( lyr[0].nBand == this.surfaceRockBand ) {
surfaceRockTemperature = lyr[0].currentTemperatureK;
// if ( ap.relativeTemperatures ) surfaceRockTemperature = lyr[0].currentTemperatureK - lyr[0].initialTemperature ;
}
// set radiative surface layer temperature
if ( lyr[0].nBand == this.thinSurfaceBand ) {
thinSurfaceBandTemperature = lyr[0].currentTemperatureK;
// if ( ap.relativeTemperatures ) thinSurfaceBandTemperature = lyr[0].currentTemperatureK - lyr[0].initialTemperature ;
// System.out.println( lyr[0].nBand + ": " + lyr[0].currentTemperatureK );
// if ( ap.year > 10000 ) lyr[0].currentTemperatureK = 322.85799999999506; // 20 deg temp jump
}
// set lowest atnospheric layer remperature
if ( lyr[0].nBand == this.airLowestBand ) {
this.airLowestBandTemperature = lyr[0].currentTemperatureK;
this.airLowestBandTemperatureK = lyr[0].currentTemperatureK;
// if ( ap.relativeTemperatures ) this.airLowestBandTemperature = lyr[0].currentTemperatureK - lyr[0].initialTemperature ;
// System.out.println( lyr[0].nBand + ": " + lyr[0].currentTemperatureK );
}
// if ( ap.icedropCounter > 0 && lyr[0].nBand == ap.icedropLayer ) {
if ( ap.icedropCounter > 0 && lyr[0].materialType == 0 ) {
// System.out.println( "Layer[" + lyr[0].nBand + "] T = " + lyr[0].currentTemperatureK );
System.out.println( "Layer[" + lyr[0].nBand + "] albedo = " + this.albedo );
ap.icedropCounter--;
}
/*
// start CO2 atmospheric warming after year 30000
if ( ap.varyAirEmissivity ) {
if ( lyr[0].genericMaterialType() == 5 && ap.year > 290000 ) {
if ( lyr[0].layerEmissivity != setAirEmissivity ) {
lyr[0].layerEmissivity = setAirEmissivity;
lyr[0].layerAbsorptivity = setAirEmissivity;
System.out.println( ap.year + ": lyr[0].layerEmissivity " + lyr[0].layerEmissivity );
}
}
}
*/
// if ( lyr[0].genericMaterialType() == 5 && ap.year < 10 ) {
// System.out.println( ap.year + ": lyr[0].layerEmissivity " + lyr[0].layerEmissivity );
// }
counter++;
}
}
nIterations = counter;
// new asteroid radius (if new layers have been added to it )
this.radius = asteroidRadius();
// report start and end of glaciation
if ( this.airMixingFactor > 0 ) {
if ( this.iceWaterDepth > 0 && this.previousIceWaterDepth == 0 ) {
// System.out.println( ap.year + " lat " + this.nLatitude + " albedo " + this.albedo + " previous albedo " + this.previousAlbedo );
this.albedoChange = true;
// change the pseudoMap latitude land to be either glaciated or unglaciated
// so that new global air temperature can be found
// System.out.println( ap.year + " glaciate" );
ap.map.glaciateLand( this.nLatitude );
} else if ( this.iceWaterDepth == 0 && this.previousIceWaterDepth > 0 ) {
this.albedoChange = true;
// System.out.println( ap.year + " deglaciate" );
ap.map.deglaciateLand( this.nLatitude );
} else {
this.albedoChange = false;
}
} else {
this.albedoChange = false;
}
this.previousIceWaterDepth = this.iceWaterDepth;
// contiguous water layers are assumed to have same temperature (due to good convective mixing)
if ( flagSetAirEmissivity && ap.year >= this.yearSetAirEmissivity ) {
setAirEmissivities( setAirEmissivity );
flagSetAirEmissivity = false;
}
if ( ap.waterRunoff ) {
// removeAllMaterialTypeLayersAboveT( 3, 273.000000001 ); // melted ice or snow
}
// removeAllMaterialTypeLayersAboveT( 3, 373.5 ); // remove steam after slight pause to rise
// remove any layer that's been flagged
if ( removeLayerFlag >= 0 ) removeLayer( removeLayerFlag );
if ( setTSLtempEqualToSRtemp ) {
setTSLtempEqualToSRtemp();
}
if ( resynchronise ) reSynchronise();
if ( flagPrintColumn ) {
printcsvColumn(0);
flagPrintColumn = false;
}
// printhcmColumn(0);
}
double surfaceTemperature( ConductiveLayer layer0, ConductiveLayer layer1, ConductiveLayer layer2, int dayIndex, int timeStepLongDays ) {
double heatflowrate = 0;
double mixedTemperature;
// radiated heat from air (and only air) in layer above
if ( layer0.genericMaterialType() == 5 ) {
// use mixedTemperature rather than layer[0].currentTemperatureK
// when airMixingFactot > 0 for the purpose of changing only heat flow from below
if ( this.airMixingFactor == 0.0 ) {
// no air mixing
heatflowrate += wattsRadiated( layer0.currentTemperatureK, layer0.layerEmissivity );
this.activeAirTemperature = layer0.currentTemperatureK;
// if ( ap.year > 9500 && ap.year < 9501 ) System.out.println( ap.year + ": " + dayIndex + " Ta " + (float)this.activeAirTemperature );
} else if ( this.airMixingFactor == 1.0 ) {
// complete air mixing
heatflowrate += wattsRadiated( ap.global.getDailyValue( dayIndex ), layer0.layerEmissivity );
this.activeAirTemperature = ap.annualGlobalMeanTemperature;
layer0.nextTemperatureK = this.activeAirTemperature; // report temperature
} else {
// partial air mixing
this.activeAirTemperature = ap.combined[ this.nLatitude ].getDailyValue( dayIndex );
if ( this.iceWaterDepth > 0 ) {
mixedTemperature = ( 1.0 - this.airMixingFactor ) * ap.combined[ this.nLatitude ].getDailyValue( dayIndex ) + ( this.airMixingFactor ) * ap.global.getDailyValue( dayIndex ) ;
// mixedTemperature = ( 1.0 - this.airMixingFactor ) * ap.snow[ this.nLatitude ].getMeanDailyValue( dayIndex-timeStepLongDays, dayIndex ) + ( this.airMixingFactor ) * ap.global.getMeanDailyValue( dayIndex-timeStepLongDays, dayIndex ) ;
if ( ap.year > 12505 && ap.year < 12506 ) System.out.println( ap.year + ": " + (float)ap.snow[ this.nLatitude ].getDailyValue( dayIndex ) + " " + (float)ap.snow[ this.nLatitude ].getMeanDailyValue( dayIndex-timeStepLongDays, dayIndex) + " ~ " + (int)this.iceWaterDepth + " " + (float)horizontalProportionOfSnow + " " + (float)ap.combined[ this.nLatitude ].getDailyValue( dayIndex ));
// if ( ap.year > 12005 && ap.year < 12006 ) System.out.println( ap.year + ": " + dayIndex + " - " + ap.global.getDailyValue( dayIndex ) + " " + ap.global.getMeanDailyValue( dayIndex-timeStepLongDays, dayIndex ) );
} else {
mixedTemperature = ( 1.0 - this.airMixingFactor ) * ap.combined[ this.nLatitude ].getDailyValue( dayIndex ) + ( this.airMixingFactor ) * ap.global.getDailyValue( dayIndex ) ;
// mixedTemperature = ( 1.0 - this.airMixingFactor ) * ap.land[ this.nLatitude ].getMeanDailyValue( dayIndex-timeStepLongDays, dayIndex ) + ( this.airMixingFactor ) * ap.global.getMeanDailyValue( dayIndex-timeStepLongDays, dayIndex ) ;
if ( ap.year > 100 && ap.year < 101 ) System.out.println( ap.year + ": " + dayIndex + " - " + ap.global.getDailyValue( dayIndex ) + " " + ap.global.getMeanDailyValue( dayIndex-timeStepLongDays, dayIndex ) );
}
heatflowrate += wattsRadiated( mixedTemperature, layer0.layerEmissivity );
this.activeAirTemperature = mixedTemperature;
layer0.nextTemperatureK = this.activeAirTemperature; // report temperature
}
// store the active air temperature used in surface temperature calculation
ap.active[ this.nLatitude ].setDailyValue( this.activeAirTemperature, dayIndex );
}
// conductive heat flow from layer below
heatflowrate += wattsConducted( layer2.currentTemperatureK, layer1.currentTemperatureK, layer2.layerConductivity, layer2.height , layer2.area );
// solar irradiation
heatflowrate += this.solarGain * ( 1.0 - this.albedo ); //
// find temperature of surface assuming emissivity of 1.0
double t = Math.pow(( heatflowrate / this.stefanBoltzmannConstant ), 0.25 );
// System.out.println( ap.year + " surface hfr3= " + (float)heatflowrate + " t " + (float)t );
return( t );
}
// air temperature over chequerboard of snow and land ranges between the two
double terrainModifiedAirTemperature( int dayIndex ) {
double ta = this.horizontalProportionOfSnow * ap.snow[ this.nLatitude ].getDailyValue( dayIndex )
+ ( 1.0 - this.horizontalProportionOfSnow ) * ap.land[ this.nLatitude ].getDailyValue( dayIndex );
return( ta );
}
// conductive heat flow
double wattsConducted( double t1, double t2, double k, double d, double a ) {
double watts = ( t1 - t2 ) * k * a / d;
return( watts );
}
// Boltzmann
double wattsRadiated( double t, double e ) {
double s = 5.6699E-8; // Stefan's constant
double watts = e * s * Math.pow( t, 4.0 );
return( watts );
}
// calculate conduction heat flow rate (Watts)
// k = conductivity, a = area, d = depth, t1-t2 temp difference
double conductiveHeatFlow( ConductiveLayer l1, double k, double a, double d, double t1, double t2 ) {
double heatflowrate = (t1 - t2) * k * a / d;
// if ( l1.nBand == 30 && l1.materialType == 3 ) System.out.println( "k a d R t1 t2 hfr " + (float)k + ", " + (float)a + ", " + (float)d + ", " + (float)(d/(k*a)) + ", " + (float)t1 + ", " + (float)t2 + ", " + (float)heatflowrate );
return ( heatflowrate );
}
// calculate conduction heat flow rate (Watts) between layer centres
// thermal resistance between layers = r1/2 + r2/2
double conductiveHeatFlow( ConductiveLayer l1, double r1, double r2, double t1, double t2 ) {
double rm = ( r1 + r2 ) / 2.0;
// System.out.println( "rm " + rm + " r1 " + r1 + " r2 " + r2 );
double heatflowrate = (t1 - t2) / rm;
// if ( l1.nBand == 30 && l1.materialType == 3 ) System.out.println( "k r1 rm t1 t2 hfr " + (float)l1.layerConductivity + ", " + (float)r1 + ", " + (float)rm + ", " + (float)t1 + ", " + (float)t2 + ", " + (float)heatflowrate );
// System.out.println( "rm " + rm + " hfr " + heatflowrate );
return ( heatflowrate );
}
// from https://en.wikipedia.org/wiki/Planck%27s_law
// freq = frequency, 1/s
double planckBlackBodyFunction( double tK, double freq ) {
double B = 0; // spectral radiance
double h = 6.626070040E-34; // Planck constant Js
double c = 299792458; // velocity of light in vacuum or air m/s
double kB = 1.38064852E-23; // Boltzmann constant J/K
double f1 = ( 2.0 * h * freq * freq * freq ) / ( c * c );
double f2 = ( h * freq ) / ( kB * tK );
B = f1 / ( Math.pow( Math.E, f2 ) - 1.0 );
return( B );
}
// source https://www.acs.org/content/acs/en/climatescience/atmosphericwarming/singlelayermodel.html
// (incoming) (1 – a)Save = (1 – e)sTp4 + esTa4 (outgoing) . . . . . . . . .(1)
// (absorbed) esTp4 = 2 esTa4 (emitted) . . . . . . . . . . . . . . . . . . (2)
// Ta4 = (1/2) Tp4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .(3)
// We can substitute for Ta4 in the planetary balance equation (1) and solve for Tp
// Tp = {[2 (1 – a)Save]/ [s(2 – e)]}1/4
// tested: results are same as in ACS
void ACSalgebra() {
double save = 342.0; // W/m2 mean solar radiation
double albedo = 0.3; // earth mean albedo
double sigma = 5.67E-8; // Stefan-Boltzmann constant, 5.67·10–8 W·m–2·K–4
double emissivity, ta;
System.out.println( "ACS single layer atmosphere" );
for ( emissivity=0; emissivity<=1.0; emissivity += 0.1 ) {
double tp = Math.pow ( ( 2.0 * (1.0 - albedo) * save ) / ( sigma * ( 2.0 - emissivity ) ), 0.25 );
ta = tp * Math.pow( 0.5, 0.25 );
System.out.println( "air emissivyty " + (float)emissivity + " surface temperature " + (float)tp + " air temperature " + (float)ta );
}
}
// set thin surface layer temperature equal to surface rock temperature in layer below
// this stops heat flow from surface rock to thin surface layer
void setTSLtempEqualToSRtemp() {
double t;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
// layer, and set fixTemperature flag
while( itr.hasNext() ){
layer = itr.next();
if ( layer.nBand == surfaceRockBand ) {
t = layer.currentTemperatureK;
// System.out.print( "layer " + layer.nBand + " temperature " + (float)t );
layer = itr.previous(); // need to go back twice to get previous layer
layer = itr.previous();
layer.currentTemperatureK = t;
// System.out.println( " becomes layer " + layer.nBand + " temperature " );
break;
}
}
}
void setInitialTemperatures() {
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
// set initial Temperature to current temperature
while( itr.hasNext() ){
layer = itr.next();
layer.initialTemperature = layer.currentTemperatureK;
}
}
void reduceLayerHeight( ConductiveLayer layer, double heatGain, double latentHeat, double layerDensity ) {
double meltMass = heatGain / latentHeat;
double meltVolume = meltMass / layerDensity;
double meltHeight = meltVolume; // in a 1m x 1m block, melt height = melt volume
layer.mass = layer.mass - meltMass;
layer.volume = layer.volume - meltVolume;
layer.height = layer.height - meltHeight;
System.out.println( "band " + layer.nBand + " layer height " + layer.height + " layer temp " + layer.currentTemperatureK + " layer init temp " + layer.initialTemperature );
if ( layer.height <=0 ) removeLayerFlag = layer.nBand;
}
// construct a growth/decay curve with time constant rc0 from ttarget temperature layer to tf final temperature
void setExponentialTemperatures( double rc0, double tf, double ttarget ) {
double ds = layerTemperatureDepth( ttarget );
double d0 = Math.log( tf - ttarget ) / tf;
System.out.println( ds + " " + d0 );
double tlayer, dl;
boolean found = false;
double depth = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
// find thin surface layer, and measure depth from there
while( itr.hasNext() ){
layer = itr.next();
if ( layer.materialType == 0 ) {
depth = 0;
} else {
depth += layer.height / 2.0;
if ( depth == ds ) {
found = true;
break;
} else {
depth += layer.height / 2.0;
}
}
}
if ( found ) {
while( itr.hasNext() ){
layer = itr.next();
depth += layer.height / 2.0;
dl = depth - ds + d0;
tlayer = exrponentialTemperature( tf, rc0, dl );
System.out.println( "depth " + depth + " temperature " + tlayer );
layer.currentTemperatureK = tlayer;
layer.nextTemperatureK = tlayer;
depth += layer.height / 2.0;
}
}
}
// function to generate a temperature approaching tmax at given depth (km)
// suggested value for rc time constant = 1000. Smaller values raise temperatures, larger lower them.
double exrponentialTemperature( double tmax, double rc, double depth ) {
return ( tmax - tmax * Math.pow ( Math.E, -depth/rc ) );
}
void fixBandTemperature( boolean fix, int nband ) {
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
// layer, and set fixTemperature flag
while( itr.hasNext() ){
layer = itr.next();
if ( layer.nBand == nband ) {
if ( fix ) {
layer.fixTemperature();
} else {
layer.unfixTemperature();
}
}
}
}
// thin surface layer albedo is usually -1 to indicate the albedo is that of the material bdenath it.
void setThinSurfaceBandAlbedo( double a ) {
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
// layer, and set fixTemperature flag
while( itr.hasNext() ){
layer = itr.next();
if ( layer.genericMaterialType() == 0 ) {
layer.albedo = a;
break;
}
}
}
// set active calculation band range true, remainder false
void setActiveCalculationLayers( int startBand, int stopBand ) {
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while( itr.hasNext() ){
layer = itr.next();
layer.setActiveCalculation( false );
if ( layer.nBand <= startBand && layer.nBand >= stopBand ) {
layer.setActiveCalculation( true );
}
}
}
// set ActiveCalculation flag in layers <= start band and >= stop band to desired state
// leaving other layers unchanged
void setActiveCalculationLayers( int startBand, int stopBand, boolean state ) {
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while( itr.hasNext() ){
layer = itr.next();
if ( layer.nBand <= startBand && layer.nBand >= stopBand ) {
layer.setActiveCalculation( state );
}
}
}
// find the depth below thin surface layer of layer at targetTemperature
double layerTemperatureDepth( double targetTemperature ) {
boolean found = false;
double depth = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
// find thin surface layer, and measure depth from there
while( itr.hasNext() ){
layer = itr.next();
if ( layer.materialType == 0 ) {
depth = 0;
} else {
depth += layer.height / 2.0;
if ( layer.currentTemperatureK >= targetTemperature ) {
break;
} else {
depth += layer.height / 2.0;
}
}
}
return( depth );
}
double volumeAtmosphere( double massAtmosphere, double densityAtmosphere ) {
double volume = 0;
volume = massAtmosphere / densityAtmosphere;
return( volume );
}
double heightAtmosphere( double volumeAtmosphere ) {
double h = 0;
double volumeColumn = this.sphereVolume( this.referenceRadius );
return( h );
}
double volumeOneSquareMetre( ) {
double volume = 0;
return( volume );
}
// area of circular strip at column latitude
double areaAtLatitude() {
double stripWidthDegrees = 90.0 / ap.nLatitudes;
double lat1 = this.columnLatitude + stripWidthDegrees / 2.0;
double lat2 = this.columnLatitude - stripWidthDegrees / 2.0;
double a =2.0 * Math.PI * this.referenceRadius * this.referenceRadius * ( Math.sin( lat1 * this.degreesToRadians ) - Math.sin( lat2 * this.degreesToRadians ) );
return( a );
}
// set all air layer emissivities and absorptivities
void setAirEmissivities( double e ) {
double e0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while( itr.hasNext() ){
layer = itr.next();
if ( layer.genericMaterialType() == 5 ) {
e0 = layer.layerEmissivity;;
layer.layerEmissivity = e;
layer.layerAbsorptivity = e;
System.out.println( ap.year + " layer " + layer.nBand + " air emissivity changed from " + e0 + " to " + e );
}
}
}
// set up a geothermal gradient from the thin surface layer down to a maximum depth
void setGeothermalGradient( double tThinSurface, double degreesPerKm, double maxDepth ) {
boolean found = false;
double depth;
double thinSurfaceRadius = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
// find thin surface layer
while( itr.hasNext() ){
layer = itr.next();
thinSurfaceRadius = layer.radius;
if ( layer.materialType == 0 ) {
layer.currentTemperatureK = tThinSurface;
layer.nextTemperatureK = tThinSurface;
found = true;
break;
}
}
if ( found ) {
while( itr.hasNext() ) {
layer = itr.next();
depth = thinSurfaceRadius - layer.radius;
depth += layer.height / 2.0; // depth of middle of layer'
depth = depth / 1000.0; // depth in km
if ( depth <= maxDepth ) {
layer.currentTemperatureK = tThinSurface + depth * degreesPerKm;
layer.nextTemperatureK = layer.currentTemperatureK;
} else {
break;
}
}
} else {
System.out.println( "thin surface layer not found " );
}
}
double terrainAlbedo( double terrainHeight, double terrainWidth, double underlyingTerrainAlbedo, double snowDepth, double snowAlbedo, int dayIndex ) {
double terrainAlbedo = underlyingTerrainAlbedo;
if ( snowDepth > 0 ) {
if ( this.iceMeltingAtBase ) {
if ( terrainHeight <= 0 ) terrainHeight = 0.001;
double snowVolume = snowDepth * terrainWidth * terrainWidth;
double valleySnowDepth = Math.sqrt( ( snowVolume * terrainHeight ) / ( terrainWidth * terrainWidth ) );
double valleySnowWidth = valleySnowDepth * terrainWidth / terrainHeight;
horizontalProportionOfSnow = valleySnowWidth / terrainWidth;
if ( horizontalProportionOfSnow > 1.0 ) horizontalProportionOfSnow = 1.0;
if ( horizontalProportionOfSnow < 0.0 ) horizontalProportionOfSnow = 0.0;
terrainAlbedo = horizontalProportionOfSnow * snowAlbedo + (1 -horizontalProportionOfSnow) * underlyingTerrainAlbedo;
// System.out.println( "horizontalProportionOfSnow " + horizontalProportionOfSnow );
} else {
// thin snow has lower albedo
if ( snowDepth > 50 ) {
terrainAlbedo = snowAlbedo;
} else {
terrainAlbedo = snowAlbedo / 2;
}
}
}
// use snow coverage to update pseudoMap and get new combined air temperature
ap.combinedAirTemperature( this.nLatitude, this.horizontalProportionOfSnow, dayIndex );
return( terrainAlbedo );
}
// s = snow depth, hf = forest height, as = albedo of snow or ice, af = albedo of forest trees
double forestAlbedo( double s, double hf, double as, double af ) {
double fAlbedo, at, ar;
// albedo of forest is albedo of underlying rock if forest height = 0
// and falls to albedo of dense forest when forest height > 5.0 m
if ( hf < 5.0 ) {
ar = 0.4; // albedo of underlying rock
at = ( 1.0 - hf / 5.0 ) * ar + ( hf / 5.0 ) * af;
} else {
at = af;
}
// at = lbedo of forest when sno height = 0
// as = albedo of snow
if ( s < hf ) {
fAlbedo = ( s / hf ) * as + ( 1.0 - s / hf ) * at;
} else {
fAlbedo = as;
}
return( fAlbedo );
}
// add up all the layer masses
double massOfAllLayers() {
double totalMass = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while( itr.hasNext() ){
layer = itr.next();
totalMass += layer.mass;
}
return( totalMass );
}
double massOfEntireAsteroid() {
double totalMass = 0;
double massSlice = this.massOfAllLayers();
double asteroidSurfaceArea = this.sphereSurfaceArea( this.referenceRadius );
totalMass = massSlice * asteroidSurfaceArea;
return( totalMass );
}
// assumes only one set of contiguous layers of water
double equaliseContiguousWaterLayerTemperatures() {
int n = 0;
double tmean;
double v = 0;
double vt = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while(itr.hasNext()){
layer = itr.next();
if ( layer.material.isWater( layer.currentTemperatureK) ) {
v += layer.volume;
vt += layer.currentTemperatureK * layer.volume;
n++;
}
}
tmean = vt / v;
if ( n > 1 ) {
itr=llayer.listIterator();
while(itr.hasNext()){
layer = itr.next();
if ( layer.material.isWater( layer.currentTemperatureK) ) {
layer.currentTemperatureK = tmean;
}
}
}
return( tmean );
}
// Measures ice presence above rock layers
// Called with control = 0, each period during which ice (mt 3 ) is present is added to total
// Called with ccontrol = 1, total/period is printed out, and total reset
// assumes constnnt period
double measureDegreeOfGlaciation( double gtotal, double period, double timestep, int control ) {
double degree = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
this.iceDepth = 0;
while( itr.hasNext() ){
layer = itr.next();
// record glaciations
if ( layer.material.isIce( layer.currentTemperatureK) ) {
degree = timestep;
iceDepth += layer.height;
}
// quit when rock or iron reached
if ( layer.materialType == 2 || layer.materialType == 1 ) {
break;
}
}
// System.out.println( degree + " " + gtotal );
gtotal += degree;
if ( control == 1 ) {
// System.out.println( "% glaciation " + ( 100.0 * gtotal / period ) + " over period "+ period );
// gtotal = 0;
}
return( gtotal );
}
// set radioactive heat produced in all layers of material type mt
void setRadioactiveHeatGeneration( double mt, double wattsPerCubicMetre ) {
double totalHeatProduced = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while( itr.hasNext() ){
layer = itr.next();
if ( layer.materialType == mt ) {
layer.heatProduced = layer.volume * wattsPerCubicMetre;
totalHeatProduced += layer.heatProduced;
}
}
System.out.println( "total radioactive heat produced " + totalHeatProduced + " W/m^3");
}
// using an array of temperatures, set layer temperatures from top down
void setLayerTemperatures( double temps[], double multiplier ) {
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
int i = 0;
while( itr.hasNext() ){
layer = itr.next();
layer.currentTemperatureK = temps[i++] * multiplier;
}
}
void restoreInitialTemperatures() {
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while( itr.hasNext() ){
layer = itr.next();
layer.restoreInitialTemperature();
}
}
// using existing temperatures, adjust layer temperatures using multiplier
void adjustLayerTemperatures( double multiplier ) {
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
int i = 0;
while( itr.hasNext() ){
layer = itr.next();
layer.currentTemperatureK = layer.currentTemperatureK * multiplier;
}
}
// recalculate all the radii of layers above startband,
// using layer volume to find new heights and radii
void recalculateRadii( int startBand ) {
int count = 0;
boolean started = false;
double subsphereRadius, subsphereVolume, subsphereArea, subsphereSliceWidth, sliceFraction, layerVolume;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
subsphereRadius = this.referenceRadius;
subsphereArea = this.sphereSurfaceArea( subsphereRadius );
subsphereVolume = this.sphereVolume( subsphereRadius );
sliceFraction = 1.0 / subsphereArea;
// System.out.println( "recalculate radii " );
// first go down the layers to the start layer
itr = llayer.listIterator();
while(itr.hasNext()){
layer = itr.next();
if ( layer.nBand == startBand ) {
started = true;
break;
}
}
if ( itr.hasPrevious() ){
layer = itr.previous();
subsphereRadius = layer.radius;
subsphereVolume = this.sphereVolume( subsphereRadius );
subsphereArea = this.sphereSurfaceArea( subsphereRadius );
subsphereSliceWidth = layer.width;
sliceFraction = subsphereSliceWidth * subsphereSliceWidth / subsphereArea;
// System.out.println( layer.nBand + "- " + layer.radius + " " + layer.width );
}
// go up the layers from the centre outwards to find the start band
while( itr.hasPrevious() ){
layer = itr.previous();
if ( started ) {
// layer mass is fixed, and layer density may be variab;e
// layer volume is calculated from layermass and density
layer.volume = layer.mass / layer.layerDensity;
layerVolume = layer.volume / sliceFraction; //
layer.radius = this.layerRadius( layerVolume, subsphereVolume );
layer.height = layer.radius - subsphereRadius;
if (layer.height < 0 ) layer.height = 0;
layer.width = layer.radius / this.referenceRadius;
layer.length = layer.radius / this.referenceRadius;
// System.out.println( layer.nBand + ": r " + layer.radius + " w " + layer.width + " v " + layer.volume);
}
subsphereRadius = layer.radius;
subsphereVolume = this.sphereVolume( subsphereRadius );
subsphereArea = this.sphereSurfaceArea( subsphereRadius );
subsphereSliceWidth = layer.width;
sliceFraction = subsphereSliceWidth * subsphereSliceWidth / subsphereArea;
}
}
// remover a single contiguous set of layers of material type1 (e.g. all ice layers in ice sheet )
// for single layer removal, set both types to be the same type.
int removeAllMaterialTypeLayersAboveT( double type1, double t ) {
boolean layerRemoved = false;
int count = 0;
int nremoved = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
if ( ap.year > 12000 ) System.out.println( ap.year + " start removeAllMaterialTypeLayersAboveT" );
double heightRemoved = 0;
itr=llayer.listIterator();
while(itr.hasNext()){
layer = itr.next();
layerRemoved = false;
// remover material in 0.1 type range (e.g. 3.00 and 3.05, ice and snow)
if ( layer.materialType >= (int)type1 && layer.materialType < ((int)type1+0.2) && layer.currentTemperatureK > t ) {
itr.remove();
heightRemoved += layer.height;
nlayers--;
nremoved++;
layerRemoved = true;
System.out.println( ap.year + " removed layer band " + (int)layer.nBand + " at " + layer.currentTemperatureK + " K " );
// start snowfall 10,000 yrs later
// ap.eventHandler.precipitstionCZYN = (int)(ap.eemianCZYN + ap.year + 5000.0);
// ap.eventHandler.precipitstionCZYN = 1000.0 * ( (int)( ap.eventHandler.precipitstionCZYN / 1000.0 ) );
System.out.println( "ap.eventHandler.precipitstionCZYN " + ap.eventHandler.precipitstionCZYN );
ap.eventHandler.insertEvent( ap.eventHandler.precipitstionCZYN );
} //
count++;
// if ( !layerRemoved && nremoved > 0 ) break;
}
// now go back up to surface, lowering radii of layers above
if ( nremoved > 0 ) {
layer = itr.previous();
while( itr.hasPrevious() ){
layer = itr.previous();
// System.out.println( "lowered layer type " + (int)layer.materialType + " by " + heightRemoved );
layer.radius -= heightRemoved;
count--;
}
}
if ( nremoved > 0 ) {
System.out.println( nremoved + " layers removed " );
// System.out.println( type1 + " number of bands " + setBandNumbers() + " nremoved " + nremoved );
reSynchronise(); // renumber layers
if ( ap.constantAddLayerDepther ) {
addLayer( ap.column[ ap.nviewLatitude ].thinSurfaceBand, 270, ap.addLayerDepth, ap.addLayerMaterial );
}
}
if ( ap.year > 12000 ) System.out.println( "end removeAllMaterialTypeLayersAboveT" );
return( nremoved );
}
// add new layer at nBand nb, of material mt, temperature t.
// adds UNDER layer nb
void addLayer( int nb, double t, double h, double mt ) {
double r, w, q;
q = 0;
ConductiveLayer layer, lyr;
ListIterator<ConductiveLayer> itr = null;
// System.out.println( "add layer" + nb );
itr=llayer.listIterator();
while(itr.hasNext()){
layer = itr.next();
r = layer.radius;
w = layer.width;
if ( layer.nBand == nb ) {
// use 5 year tick period
lyr = new ConductiveLayer( this, r+h, w, h, t, mt, 1.0, 0 );
itr.add( lyr );
System.out.println( ap.year + " added layer of material " + mt + " at band" + nb );
System.out.println( "layer temperature " + lyr.currentTemperatureK + " heat content per kg" + lyr.layerHeatContentPerKg );
// q = lyr.pt.temperature( lyr.layerHeatContentPerKg );
// System.out.println( "layer temperature from heat content " + q );
// System.out.println( "layer heat content from temperature " + lyr.pt.heatContentPerKg( q ) );
}
}
reSynchronise();
// printhcmColumn( 0 );
}
// Add volume v at temperature t into nband nb
// All higher layers need their height recalculated (but havent yet)
void addIntopLayer( double t, double v, int nb ) {
double r, w, vnew, tnew;
ConductiveLayer layer, lyr;
ListIterator<ConductiveLayer> itr = null;
itr = llayer.listIterator();
if (itr.hasNext() ) {
layer = itr.next();
if ( layer.nBand == nb ) {
vnew = layer.volume + v;
tnew = meanTemperature( layer.volume, layer.currentTemperatureK, v, t );
layer.volume = vnew;
layer.currentTemperatureK = tnew;
layer.height = layer.volume / ( layer.width * layer.width ); // this needs to be correctly calculated
}
}
reSynchronise();
}
// look for H2O layer (materialType 3.xxx)
// and if found add more ice/snow/water
// else create new layer
void addIceLayer( double t, double h, double mt ) {
double r, w, h1, h2, m1, m2;
ConductiveLayer layer, lyr;
ListIterator<ConductiveLayer> itr = null;
// System.out.println( "start addIceLayer" );
itr=llayer.listIterator();
while(itr.hasNext()){
layer = itr.next();
r = layer.radius;
w = layer.width;
if ( iceWaterDepth == 0 ) {
// if no water layer, create one beneath thin surface layer
if ( layer.materialType == 0 ) {
lyr = new ConductiveLayer( this, r+h, w, h, t, mt, 1.0, 0 );
itr.add( lyr );
// System.out.println( "lyr.layerHeatContent " + lyr.layerHeatContent );
ap.addedIceMass = lyr.mass; // mass to be added to grow layer
ap.addedIceHeatContentPerKg = lyr.layerHeatContentPerKg; // heat content per Kg of added snow
layer.layerHeatAdded = lyr.layerHeatContent; // added heat isn't just heat flow
// System.out.println( ap.year + " added brand new layer of snow/ice " + mt + " at band" + layer.nBand );
// System.out.println( "brand new addedIceMass " + ap.addedIceMass + " addedIceHeatContent " + ap.addedIceHeatContent );
}
} else {
// if water layer, add to it, finding new mean temperature
// assumes adding H2O to H2O
if ( layer.materialType >= 3.0 && layer.materialType < 3.2 ) {
// find new layer mass and temperature
// layer.mass += ap.addedIceMass;
// layer.layerHeatAdded += ap.addedIceHeatContent; // added heat isn't just heat flow
// layer.layerHeatContentPerKg += ap.addedIceHeatContent / layer.mass;
// layer.layerHeatContent = layer.layerHeatContentPerKg * layer.mass;
m1 = layer.mass;
m2 = ap.addedIceMass;
h1 = layer.layerHeatContentPerKg;
h2 = ap.addedIceHeatContentPerKg;
layer.layerHeatContentPerKg = ( h1 * m1 + h2 * m2 ) / ( m1 + m2 );
layer.mass = ( m1 + m2 );
layer.layerHeatContent = layer.layerHeatContentPerKg * layer.mass;
// System.out.println( "layer.layerHeatContent" + (float)layer.layerHeatContent + " added mass " + (float)ap.addedIceMass + " heatcontent " + (float)ap.addedIceHeatContent );
layer.currentTemperatureK = layer.pt.temperature( layer.layerHeatContentPerKg );
// redimension layer and layer above
// this.recalculateRadii( layer.nBand );
//
// System.out.println( ap.year + " new snow/ice layer heat content " + (float)layer.layerHeatContent + " temperature " + (float)layer.currentTemperatureK);
break;
}
}
}
// System.out.println( "resynchronise" );
reSynchronise();
// printhcmColumn( 71 );
// System.out.println( "end addIceLayer" );
}
// remove nBand layer nb;
void removeLayer( int nb ) {
double snowMass = 0;
double pe = 0;
int uni = 0;
int nremoved = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while( itr.hasNext() ){
layer = itr.next();
if ( layer.genericMaterialType() == 3 ) snowMass += layer.mass;
if ( layer.nBand == nb ) {
if ( layer.genericMaterialType() == 3 ) {
// add potential energy of overlying snow to layer below removed snow layer
pe = snowMass * layer.height * 9.81;
if ( pe > 0 ) addHeatToLayer( nb-1, pe );
// layerAbove.layerHeatContent += pe;
// layerAbove.layerHeatContentPerKg = layerAbove.layerHeatContent / layerAbove.mass;
// System.out.println( "P.E. added " + (float)pe + " layer heat content " + (float)layerAbove.layerHeatContent );
}
itr.remove();
nlayers--;
nremoved++;
}
// layerAbove = layer;
}
// System.out.println( ap.year + " band " + nb + ": " + uni + " removed on request ");
removeLayerFlag = -1; // flag no layer to be removed
resynchronise = true;
}
// add heat to layer nb
void addHeatToLayer( int nb, double heat ) {
ConductiveLayer layer, layerAbove;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
layerAbove = itr.next();
while( itr.hasNext() ){
layer = itr.next();
if ( layer.nBand == nb ) {
// add potential energy of overlying snow to layer bn
layer.layerHeatContent += heat;
layer.layerHeatContentPerKg = layer.layerHeatContent / layer.mass;
// System.out.println( "P.E. added " + (float)heat + " layer heat content " + (float)layer.layerHeatContent );
}
}
}
// find band number of layer of higher densiy than material mt
int findHigherDensityNBand( double mt ) {
int nb = -1;
Material mtrl = new Material( this, mt );
ConductiveLayer layer, lyr;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while(itr.hasNext()){
layer = itr.next();
System.out.println( layer.nBand + ": " + layer.layerDensity + " > " + mtrl.density + " " + nb );
if ( layer.layerDensity >= mtrl.density ) {
nb = layer.nBand;
break;
}
}
return( nb );
}
double meanTemperature( double v1, double t1, double v2, double t2 ) {
double tmean = ( v1 * t1 + v2 * t2 ) / ( v1 + v2 );
return( tmean );
}
// find mean temperature of whole asteroid
// this is the temperature that surface rocks would reach beneath a perfect insulator
double meanAsteroidTemperature() {
double vtot = 0;
double vttot = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while(itr.hasNext()){
layer = itr.next();
vttot += layer.volume * layer.currentTemperatureK;
vtot += layer.volume;
}
vttot = vttot / vtot;
return( vttot );
}
double radiusConstantVolumeLayer( int layer, double layerVolume ) {
double lRadius, vol;
double f = 4.0 * Math.PI / 3;
double cubeRoot = 1.0 / 3.0;
vol = (layer+1) * layerVolume;
lRadius = Math.pow( ( vol / f ), cubeRoot );
return( lRadius);
}
// return the radius of the innermost layer of a sphere of radius r with n layers in it.
double equalLayerVolumeSphere( int numLayers, double sphereRadius ) {
double r1;
double r[] = new double[ numLayers ];
// In a sphere with centre layer radius R1,
// the radius of the nth layer, Rn, is given by Rn^3 = n * R1^3
r1 = Math.pow( ( sphereRadius * sphereRadius * sphereRadius / numLayers ), 0.333333 );
return( r1 );
}
double sphereVolume( double r ) {
double volume = 1.3333333333 * Math.PI * r * r * r;
return( volume );
}
double sphereSurfaceArea( double r ) {
double a = 4.0 * Math.PI * r * r;
return( a );
}
double sphereRadius( double vol ) {
// sphere volume = 4/3 .pi. r^3
double fn = ( vol * 3.0 ) / ( 4.0 * Math.PI );
double r = Math.pow( fn, 0.33333333 );
return( r );
}
// find outer radius of layer sitting above an inner subsphere
double layerRadius( double layerVolume, double volumeSubsphere ) {
double radius = 0;
double fn = 3.0 * ( layerVolume + volumeSubsphere ) / ( 4.0 * Math.PI );
radius = Math.pow( fn, 1.0/3.00 );
return( radius );
}
// find highest band number of nmaterial layers
int countMaterialBands( int nmaterial ) {
int count = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while(itr.hasNext()){
layer = itr.next();
if ( layer.materialType == nmaterial ) {
count++;
}
}
return( count );
}
double findIcePlusWaterDepth() {
double depth = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while(itr.hasNext()){
layer = itr.next();
if ( layer.material.isIce( layer.currentTemperatureK) || layer.material.isWater( layer.currentTemperatureK) ) {
// System.out.println( layer.nBand + " " + layer.materialType );
depth += layer.height;
}
}
// System.out.println( "ice depth " + depth );
return( depth );
}
// find highest band number of nmaterial layers
int materialBandNumber( int nmaterial ) {
int count = -1;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while(itr.hasNext()){
layer = itr.next();
if ( layer.materialType == nmaterial ) {
count = layer.nBand;
break;
}
}
return( count );
}
// layers are counted from top surface (layer 0) of asteroid downwards.
// bands are counted from asteroid centre (band 0) upwards.
// This needs to be called whenever layers are added to or subtracted from the asteroid
// This should also be called fairly regularly to update graphic display temperatures and materials in layers
int setBandNumbers() {
boolean flag1 = false;
boolean flag2 = false;
boolean flag3 = false;
int n = 0;
int count = nlayers - 1;
double totalMass = 0;
this.topLayerOfAtmosphere = count;
iceWaterDepth = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
// seems that list iterator can only go back to revious elemnts if it's been forrward to next elemnts.
while( itr.hasNext() ){
layer = itr.next();
// store data for graphics display
graphicLayer[n][0] = layer.nBand;
graphicLayer[n][1] = layer.materialType;
graphicLayer[n][2] = layer.radius;
graphicLayer[n][3] = layer.height;
graphicLayer[n][4] = layer.currentTemperatureK;
graphicLayer[n][5] = layer.phase;
graphicLayer[n][6] = layer.upwardHeatFlowRate;
graphicLayer[n][7] = layer.uniqueNumber;
graphicLayer[n][8] = layer.material.colourIndex[ layer.phase ];
graphicLayer[n][9] = layer.initialTemperature;
layer.nBand = count;
totalMass += layer.mass;
// find thin surface band for solar irradiation
if ( layer.materialType == 0 && !flag1 ) {
this.thinSurfaceBand = layer.nBand;
this.surfaceBand = layer.nBand - 1; // old surface band is band underneath thin surface band
// System.out.println( "thin layer found " + layer.nBand );
flag1 = true;
}
// find top rock layer
if ( layer.materialType == 2 && !flag2 ) {
this.surfaceRockBand = layer.nBand;
// System.out.println( "surf layer found " + layer.nBand );
flag2 = true;
}
// record temperature of surface band under thin surface band
if ( layer.nBand == surfaceBand ) {
this.surfaceBandTemperature = layer.currentTemperatureK;
}
// find ice + water depth
if ( layer.materialType >= 3.0 && layer.materialType < 3.2 ) {
iceWaterDepth += layer.height;
}
if ( layer.genericMaterialType() == 5 && !flag3 ) {
this.airLowestBand = layer.nBand;
// System.out.println( "sair bottom band " + layer.nBand );
// flag3 = true;
}
count--;
n++;
}
// set topmost layer radius to pre-defined top of atmospher
// graphicLayer[n-1][2] = this.topOfAtmosphere;
graphicLayerCount = n;
graphicLayer[n][0] = -1; // terminator
layer = itr.previous();
graphicSurfaceRockBand = this.surfaceRockBand;
// System.out.println( "band thinSurf " + this.thinSurfaceBand + " surf " + this.surfaceBand + " rock " + this.surfaceRockBand );
// System.out.println( "setBandNumbers totalMass " + totalMass );
return( n );
}
// add up layer heights to get asteroid radius
// and set nlayers variable
double asteroidRadius() {
int count = 0;
double r = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while(itr.hasNext()){
layer = itr.next();
r += layer.height;
count++;
}
this.nlayers = count;
return( r );
}
// adapted from GlobalGlaciation
double AreaAtLatitude( double lat, double latRange ) {
double a;
double l1 = degreesToRadians * ( lat + latRange );
double l2 = degreesToRadians * ( lat - latRange );
if ( lat < 90 ) {
a = 2.0 * Math.PI * referenceRadius * referenceRadius * ( Math.sin( l1 ) - Math.sin( l2 ) );
} else {
a = 2.0 * Math.PI * referenceRadius * referenceRadius * ( 1.0 - Math.sin( degreesToRadians * (90.0 - latRange/2.0) ) );
// System.out.println( Math.sin( degreesToRadians * latRange/2.0 ) );
// a = 0;
}
// System.out.println( lat + " " + a );
return( a );
}
void setLayerTicks( int startBand, int count1 ) {
int count = 0;
double r = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while(itr.hasNext()){
layer = itr.next();
if ( layer.nBand < startBand ) {
count = layer.setupTick( count1 );
System.out.println( layer.nBand + " tick count " + count );
}
count++;
}
}
// find layer pressures, and use these and layer temperatures to find layer strain,
// and then use layer strain to find layer density and volume
// and use layer volume to find layer width, height, and radius from centre
// ( used to only calculate rock pressures. It doesn't redimension layers.)
void redimensionLayers() {
double glocal; // local gravitational acceleration
double layerCentreRadius;
double totalMass = 0;
double totalForce = 0;
double midLayerForce, deltaT, deltaP, depth;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while ( itr.hasNext() ) {
layer = itr.next();
layerCentreRadius = layer.radius - layer.height / 2.0;
depth = this.nominalRadius - layer.radius;
if ( depth >= 0 ) {
glocal = gravitationalAccelerationInsideAsteroid( layerCentreRadius, this.nominalRadius, this.surfaceGravity );
} else {
glocal = surfaceGravity; // atmosphere g
}
midLayerForce = totalForce + glocal * layer.mass / 2.0;
// System.out.println( layer.nBand + " midLayerForce " + midLayerForce + " totalForce " + totalForce + " glocal * layer.mass / 2.0 " + glocal * layer.mass / 2.0);
layer.previousLayerPressure = layer.currentLayerPressure;
layer.currentLayerPressure = midLayerForce / layer.area;
deltaP = layer.currentLayerPressure - layer.previousLayerPressure;
deltaT = layer.currentTemperatureK - layer.previousTemperaureK;
totalForce += layer.mass * glocal;
totalMass += layer.mass;
System.out.println( layer.nBand + " area " + layer.area + " km mass "+ (int)layer.mass + " kg pressure " + (float)(layer.currentLayerPressure/1000000000.0) + " gPascals" + " glocal " + (float)glocal );
}
System.out.println( "totalMass " + (float)totalMass * this.surfaceArea + " kg" );
}
void findNonTaperingLayerPressures( double a ) {
double glocal; // local gravitational acceleration
double layerCentreRadius;
double totalMass = 0;
double totalForce = 0;
double midLayerForce, deltaT, deltaP, depth, mass;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while ( itr.hasNext() ) {
layer = itr.next();
layerCentreRadius = layer.radius - layer.height / 2.0;
depth = this.nominalRadius - layer.radius;
glocal = gravitationalAccelerationInsideAsteroid( layerCentreRadius, this.nominalRadius, this.surfaceGravity );
mass = a * layer.height * layer.layerDensity;
midLayerForce = totalForce + glocal * mass / 2.0;
// System.out.println( layer.nBand + " midLayerForce " + midLayerForce + " totalForce " + totalForce + " glocal * layer.mass / 2.0 " + glocal * mass / 2.0);
layer.previousLayerPressure = layer.currentLayerPressure;
layer.currentLayerPressure = midLayerForce / a ;
deltaP = layer.currentLayerPressure - layer.previousLayerPressure;
deltaT = layer.currentTemperatureK - layer.previousTemperaureK;
totalForce += mass * glocal;
totalMass += mass;
System.out.println( layer.nBand + " km mass "+ mass + " kg pressure " + (layer.currentLayerPressure/1000000000.0) + " gPascals" + " glocal " + (float)glocal );
}
System.out.println( "totalMass " + (float)totalMass * this.surfaceArea + " kg" );
}
// layer strain (change of height) due to both pressure and temperature
double layerStrain() {
double deltaH = 0;
return( deltaH );
}
// source http://physicsteacher.in/2017/10/18/acceleration-due-to-gravity-height-depth/
// aradius = asteroid radius, r = required radius
double gravitationalAccelerationInsideAsteroid( double r, double aradius, double surfaceacceleration ) {
// System.out.println( r + ", " + aradius + ", " + surfaceacceleration );
double h = aradius - r;
return( surfaceacceleration * ( 1.0 - h / aradius ) );
}
void findLayerPressures() {
double glocal; // local gravitational acceleration
double layerCentreRadius;
double totalMass = 0;
double totalForce = 0;
double midLayerForce, deltaT, deltaP, depth;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
while ( itr.hasNext() ) {
layer = itr.next();
layerCentreRadius = layer.radius - layer.height / 2.0;
depth = this.nominalRadius - layer.radius;
if ( depth >= 0 ) {
glocal = gravitationalAccelerationInsideAsteroid( layerCentreRadius, this.nominalRadius, this.surfaceGravity );
} else {
glocal = surfaceGravity; // atmosphere g
}
midLayerForce = totalForce + glocal * layer.mass / 2.0;
// System.out.println( layer.nBand + " midLayerForce " + midLayerForce + " totalForce " + totalForce + " glocal * layer.mass / 2.0 " + glocal * layer.mass / 2.0);
layer.previousLayerPressure = layer.currentLayerPressure;
layer.currentLayerPressure = midLayerForce / layer.area;
deltaP = layer.currentLayerPressure - layer.previousLayerPressure;
deltaT = layer.currentTemperatureK - layer.previousTemperaureK;
totalForce += layer.mass * glocal;
totalMass += layer.mass;
// System.out.println( layer.nBand + " area " + layer.area + " km mass "+ (int)layer.mass + " kg pressure " + (float)(layer.currentLayerPressure/1000000000.0) + " gPascals" + " glocal " + (float)glocal );
}
// System.out.println( "totalMass " + (float)totalMass * this.surfaceArea + " kg" );
}
// when numbers of layers change, or height/thickness of layers changes,
// all the layer radii and band numbers must be recalcu;ated
void reSynchronise() {
// System.out.println( "resynchronise "+ (this.ap.currentCZYN - this.ap.eemianCZYN) );
this.radius = asteroidRadius(); // get complete radius of asteroid using all layers
this.topOfAtmosphere = this.radius;
setBandNumbers();
findLayerPressures();
recalculateRadii( surfaceRockBand );
// radius of each layer is radius of top of layer
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr = llayer.listIterator();
// double heightOfLayersAbove = 0;
while ( itr.hasNext() ){
layer = itr.next();
// layer.radius = this.radius - heightOfLayersAbove;
// heightOfLayersAbove += layer.height;
this.layerTable[ layer.nBand ] = layer;
// if ( layer.nBand > 77 ) System.out.println( ap.year + " this.layerTable[" +layer.nBand+ "].radius " + this.layerTable[ layer.nBand ].radius );
}
// if ( this.materialBandNumber( 4 ) > this.surfaceBand ) this.surfaceBand = this.materialBandNumber( 4 );
// if ( this.materialBandNumber( 3 ) > this.surfaceBand ) this.surfaceBand = this.materialBandNumber( 3 );
resynchronise = false; // unflag reSynchroise() call needed
}
void printWholeColumn() {
double mass[] = new double[10];
String name[] = new String[10];
int index;
double wLayer, r = 0;
double totalMass = 0;
double hfl = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
System.out.println( "Year " + ap.year);
// seems that list iterator can only go back to revious elemnts if it's been forrward to next elemnts.
while( itr.hasNext() ){
layer = itr.next();
index = (int)layer.materialType;
name[ index ] = layer.material.name;
mass[ index ] += layer.mass;
}
for( index=0; index<=9; index++ ) {
System.out.println( index + " " + name[index] + " mass " + mass[index] * this.surfaceArea );
}
}
void printColumn( int ref ) {
double wLayer, r = 0;
double totalMass = 0;
double totalCapacity = 0;
double airMass = 0;
double hfl = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
System.out.println( "Year " + ap.year);
// seems that list iterator can only go back to revious elemnts if it's been forrward to next elemnts.
while( itr.hasNext() ){
layer = itr.next();
totalMass += layer.mass;
totalCapacity += layer.mass * layer.layerSpecificHeat;
wLayer = layer.radius / referenceRadius;
if ( layer.genericMaterialType() == 5 ) airMass += layer.mass;
if ( ref == 1 || ref == 2 ) layer.storedTemperature[ ref-1 ] = layer.currentTemperatureK;
System.out.println( "band " + layer.nBand + " r " + (float)layer.radius + " ht " + (float)layer.height + " lvol " + (float)layer.volume + " mtl " + (float)layer.materialType + " mass " + (int)layer.mass + " lres " + (float)layer.layerResistance + " uhfr " + (float)layer.upwardHeatFlowRate + " degrees K " + (float)layer.currentTemperatureK ) ;
// System.out.println( "band " + layer.nBand + " radius " + (float)layer.radius + " height " + (float)layer.height + " lvolume " + (float)layer.layerVolume( layer.radius, layer.radius - layer.height, wLayer ) + " material " + (float)layer.materialType + " heat " + (float)layer.heatProduced + " degrees K " + (float)layer.currentTemperatureK ) ;
r += layer.height;
hfl += layer.upwardHeatFlowRate;
}
System.out.println( "asteroid radius " + r + " net upward heat flow rate " + hfl );
System.out.println( "asteroid column mass " + totalMass + " asteroid mass " + totalMass * this.surfaceArea );
System.out.println( "True Earth mass = 5.972 × 10^24 kg");
System.out.println( "Asteroid atmosphere mass = " + (float)(airMass * this.surfaceArea) + " true value 5.14E18 Kg");
System.out.println( "asteroid heat content " + this.asteroidHeatContent + " J");
System.out.println( "asteroid mean temperature " + this.asteroidHeatContent/totalCapacity + " deg K ");
}
// print .csv layer data
void printcsvColumn( int ref ) {
double heatContent = 0;
double v1 = 0;
double v2 = 0;
double v3 = 0;
double v4 = 0;
double v5 = 0;
double v6 = 0;
double wLayer, r = 0;
double totalMass = 0;
double totalCapacity = 0;
double hfl = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
Date date = new Date();
SimpleDateFormat sdf = new SimpleDateFormat("HH:mm:ss.SSS");
String str = sdf.format(date);
System.out.println( "Year " + ap.year + " " + str + " latitude " + this.columnLatitude + " " + this.test );
// seems that list iterator can only go back to revious elemnts if it's been forrward to next elemnts.
while( itr.hasNext() ){
layer = itr.next();
if ( layer.materialType == 3.05 ) {
heatContent = layer.layerHeatContent;
v1 = layer.layerConductivity;
v2 = layer.layerResistance;
v3 = layer.material.conductivity;
v4 = layer.layerSpecificHeat;
v5 = layer.material.specificHeat;
v6 = layer.pt.FractionPhaseChange;
}
totalMass += layer.mass;
totalCapacity += layer.mass * layer.layerSpecificHeat;
wLayer = layer.radius / referenceRadius;
if ( ref == 1 || ref == 2 ) layer.storedTemperature[ ref-1 ] = layer.currentTemperatureK;
System.out.println( layer.nBand + ", " + layer.materialType + ", " + (float)layer.radius + ", " + (float)layer.height + ", " + (float)layer.width + ", " + (float)layer.upwardHeatFlowRate + ", " + (float)layer.currentTemperatureK + ", " ) ;
// System.out.println( layer.nBand + ", " + layer.uniqueNumber + ", " + layer.materialType + ", " + layer.radius + ", " + layer.height + ", " + layer.width + ", " + layer.pseudoTemperatureK + ", " + layer.currentTemperatureK + ", " ) ;
r += layer.height;
hfl += layer.upwardHeatFlowRate;
}
System.out.println( "asteroid radius " + r + " net upward heat flow rate " + hfl );
System.out.println( "asteroid column mass " + totalMass + " asteroid mass " + totalMass * this.surfaceArea );
System.out.println( "True Earth mass = 5.972 × 10^24 kg");
System.out.println( "asteroid heat content " + this.asteroidHeatContent + " J");
System.out.println( "asteroid mean temperature " + this.asteroidHeatContent/totalCapacity + " deg K ");
System.out.println( "Snow layer heat content " + heatContent + " " + + v1 + " " + v2 + " " + v3 + " " + v4 + " " + v5 + " " + v6);
double TotalTime = System.currentTimeMillis() - ap.stime;
System.out.println( "Elapsed milliseconds " + TotalTime );
System.out.println( (ap.currentCZYN - ap.eemianCZYN) + " " + (ap.eventHandler.precipitstionCZYN - ap.eemianCZYN) );
System.out.println( "H20LayerHeatContent " + ap.H20LayerHeatContent + " " + ap.H20LayerMassxSpht + " " + ap.H20LayerMassxLatht + " added " + ap.H20LayerHeatAdded );
}
// print .csv layer data (re-copied from Glaciation140)
void printhcmColumn( int ref ) {
double heatContent = 0;
double wLayer, r = 0;
double totalMass = 0;
double totalCapacity = 0;
double hfl = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
Date date = new Date();
SimpleDateFormat sdf = new SimpleDateFormat("HH:mm:ss.SSS");
String str = sdf.format(date);
System.out.println( "Year " + ap.year + " " + str + " latitude " + this.columnLatitude + " " + " .hcm " );
System.out.println( "Yr " + " mtrl " + " radius " + " height " + " mass " + " timestep " + " heat/kg " + " temp K " );
// seems that list iterator can only go back to revious elemnts if it's been forrward to next elemnts.
while( itr.hasNext() ){
layer = itr.next();
totalMass += layer.mass;
totalCapacity += layer.mass * layer.layerSpecificHeat;
wLayer = layer.radius / referenceRadius;
System.out.println( layer.nBand + ", " + layer.materialType + ", " + (float)layer.radius + ", " + (float)layer.height + ", " + (float)layer.mass + ", " + (float)layer.tickDuration + ", " + (float)layer.layerHeatContentPerKg + ", " + (float)layer.currentTemperatureK + ", " ) ;
r += layer.height;
// hfl += layer.upwardHeatFlowRate;
if ( layer.nBand == ref ) break;
}
System.out.println( "thinSurfaceBand " + this.thinSurfaceBand + " thinSurfaceLayerAlbedo " + thinSurfaceLayerAlbedo );
}
// print asteroid R, C, and RC time constants
void printColumnRC( int ref ) {
double wLayer, r = 0;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr=llayer.listIterator();
System.out.println( "Year " + ap.year);
// seems that list iterator can only go back to revious elemnts if it's been forrward to next elemnts.
while( itr.hasNext() ){
layer = itr.next();
wLayer = layer.radius / referenceRadius;
if ( ref == 1 || ref == 2 ) layer.storedTemperature[ ref-1 ] = layer.currentTemperatureK;
System.out.println( "band " + layer.nBand + " v " + (float)layer.volume + " mtl " + (float)layer.materialType + " l_k " + (float)layer.layerConductivity + " lres " + (float)layer.layerResistance + " cap " + (float)layer.layerThermalCapacity + " time const (years) " + (float)(layer.layerTimeConstant / ap.secondsPerYear) ) ;
// System.out.println( "band " + layer.nBand + " radius " + (float)layer.radius + " height " + (float)layer.height + " lvolume " + (float)layer.layerVolume( layer.radius, layer.radius - layer.height, wLayer ) + " material " + (float)layer.materialType + " heat " + (float)layer.heatProduced + " degrees K " + (float)layer.currentTemperatureK ) ;
r += layer.height;
}
System.out.println( "asteroid radius " + r );
}
void printGraphicLayers() {
System.out.println( "Year " + ap.year + " graphic layers");
for ( int n=0; n<graphicLayerCount; n++ ) {
System.out.println( "band " + (int)graphicLayer[n][0] + " type " + (int)graphicLayer[n][1] + " radius " + graphicLayer[n][2] );
}
}
}
class ConductiveLayer {
GeologicalColumn astro;
PhaseTemperature pt;
double materialType; // 1=iron, 2= rock, 3=ice.water, 5=air, 6.nnn=rock-ice where .nnn = fgarnite
Material material;
double radius; // radial distance of top of layer from asteroid centre
double length; // longitudinal m
double width; // latitudinal m
double height; // radial thickness m
double area; // length x widthe
double volume; // length x width * height
double mass; // layer mass kg
double albedo; // layer albedo is set to material albedo, but can change
int phase; // 0 solid, 1 liqid. 2 vapour
double volumetricFractionGranite;
double heatProduced; // watts heat generated by radioactive granite in layer
double upwardHeatFlowRate; // conductive heat flow rate to layer above
double maxUpwardHeatFlowRate = 0;
double netHeatFlowRate; // layer heat flow from all sources
double previousTemperaureK; // temperature at previous iteration
double currentTemperatureK;
double initialTemperature; // temperature at initiation
double nextTemperatureK;
double nextPhase;
boolean fixTemperature = false; // allow or disallow temperature change
double pseudoTemperatureK; // temperature used during phase change
double phaseTemperatureK; // temperature of current phase change
boolean phaseChange; // flagged true if currently changing phase
boolean phaseHasChanged = false; // flag set when material phase changes
boolean layerPhase; // 0 solide, 1 liqid, 2 vapour, -1 phasechange
double latentHeat;
double layerDensity;
double layerConductivity;
double layerSpecificHeat;
double layerLatentHeatOfFusion;
double layerEmissivity;
double layerAbsorptivity;
double layerHeatContent; // J
double layerHeatContentPerKg; // J / Kg
double initialLayerHeatContent;
double initialLayerHeatContentPerKg;
double layerHeatAdded; // summed heat gains
double layerThermalCapacity;
double layerResistance;
double layerTimeConstant;
double currentLayerPressure; // kg / square metre
double previousLayerPressure; // kg / square metre
double layerTimestep; // layer timestep (seconds)
double tickDuration; // duration of ap tick (seconds)
int tickCounter; // current layer tick countdown
int tickResetCount; // reset value for tickCounter;
int nBand; // band numbers have reverse sense of layers, Band 0 is at centre.
double storedTemperature[] = new double[2];
double phaseSum = 0;
boolean updated;
int uniqueNumber = 0;
boolean activeCalculation = true;
Average averageTemperatureK;
Material material0;
Material material1;
Material material2;
double spHt0;
double spHt1;
double spHt2;
// default constructor
ConductiveLayer() { }
// indexed material table driven constructor (called by csvreader()) Sep 2019
ConductiveLayer( GeologicalColumn a, double rlyr, double w, double h, double t, double mt, double tickYears, double wattsPerCubicMetre ) {
this.astro = a;
getUniqueLayerNumber();
this.materialType = mt;
this.length = w;
this.width = w;
this.height = h;
this.area = w * w;
this.volume = this.layerVolume( rlyr, rlyr - this.height, this.width );
this.radius = rlyr; // radial distance of top of layer from asteroid centre
this.currentTemperatureK = t;
this.nextTemperatureK = t;
this.initialTemperature = t;
this.material = new Material( astro, mt );
// this.heatProductionGranite = 0; // watts / kg (multiple of 0.0000000001 )
this.albedo = material.albedo;
this.layerDensity = material.density;
this.layerConductivity = material.conductivity;
this.layerSpecificHeat = material.specificHeat;
this.layerLatentHeatOfFusion = material.latentHeatOfFusion;
this.layerEmissivity = material.emissivity;
this.layerAbsorptivity = material.absorptivity;
this.mass = this.volume * this.layerDensity;
this.layerResistance = this.layerResistance();
this.layerThermalCapacity = this.layerThermalCapacity();
this.layerTimeConstant = this.layerResistance * this.layerThermalCapacity;
this.heatProduced = volume * wattsPerCubicMetre;
this.setLayerPhase();
this.setTickYears( tickYears );
// System.out.println( "tickYears " + tickYears + " tick duration " +this.tickDuration + " tick count " + this.tickResetCount );
// System.out.println( "Layer constructor current temperature " + this.currentTemperatureK + " heat content " + this.layerHeatContent );
// System.out.println( "Layer constructor tick duration " + this.tickDuration / astro.secondsPerYear + " count " + this.tickResetCount );
// phase change data
double mt0 = (int)this.materialType; // solid (ice)
double mt1 = mt0 + 0.1; // liquid (water)
double mt2 = mt0 + 0.2; // gas (water vapour)
this.material0 = new Material( astro, mt0 );
this.material1 = new Material( astro, mt1 );
this.material2 = new Material( astro, mt2 );
// double spHt0 = this.layerSpecificHeat; // use current specific heat
this.spHt0 = material0.specificHeat;
this.spHt1 = material1.specificHeat;
this.spHt2 = material2.specificHeat;
pt = new PhaseTemperature( this );
System.out.println( " initial temperature " + pt.temperature( layerHeatContentPerKg ) );
this.layerHeatContentPerKg = pt.heatContentPerKg( this.currentTemperatureK );
this.layerHeatContent = this.layerHeatContentPerKg * this.mass;
this.initialLayerHeatContent = this.layerHeatContent;
this.initialLayerHeatContentPerKg = this.layerHeatContentPerKg;
this.layerHeatAdded = 0;
averageTemperatureK = new Average( (int)( 1.0 / astro.ap.timeStepYears ) );
// if ( mt == 0 ) System.out.println( " init m " + (float)this.mass + ", dy " + (float)this.layerDensity + ", v " + (float)this.volume + ", h " + (float)this.height + ", w " + (float)this.width + ", a " + (float)this.area + ", r " + (float)this.radius );
}
// indexed material table driven constructor Mar 2019 (called by hcmreader())
// System.out.println( layer.nBand + ", " + layer.materialType + ", " + (float)layer.radius + ", " + (float)layer.height + ", " + (float)layer.mass + ", " + (float)layer.layerHeatContentPerKg + ", " + (float)layer.currentTemperatureK + ", " ) ;
// band numbers, radii, height, and temperature are all recalculated from mass, density, and heat content (and so are simply included )
// band number material radius height mass heat content temperature
ConductiveLayer( GeologicalColumn a, int band, double mt, double r, double h, double m, double tick, double heatPerKg, double t ) {
double tickYears = 0; // time step one day
this.astro = a;
getUniqueLayerNumber();
this.materialType = mt;
// claculate layer dimenssions from layer mass and material density
this.material = new Material( astro, mt );
this.layerDensity = material.density;
this.mass = m;
this.volume = layerVolume( this.mass, this.layerDensity );
this.radius = radiusTopOfLayer( this.volume, r, a.referenceRadius );
setLayerDimensions( this.radius, r, a.referenceRadius );
// if ( mt == 0 ) System.out.println( " init m " + (float)this.mass + ", dy " + (float)this.layerDensity + ", v " + (float)this.volume + ", h " + (float)this.height + ", w " + (float)this.width + ", a " + (float)this.area + ", r " + (float)this.radius );
// this.heatProductionGranite = 0; // watts / kg (multiple of 0.0000000001 )
this.albedo = material.albedo;
this.layerConductivity = material.conductivity;
this.layerSpecificHeat = material.specificHeat;
this.layerLatentHeatOfFusion = material.latentHeatOfFusion;
this.layerEmissivity = material.emissivity;
this.layerAbsorptivity = material.absorptivity;
this.mass = this.volume * this.layerDensity;
this.layerResistance = this.layerResistance();
this.layerThermalCapacity = this.layerThermalCapacity();
this.layerTimeConstant = this.layerResistance * this.layerThermalCapacity;
// this.heatProduced = volume * wattsPerCubicMetre;
this.layerHeatContentPerKg = heatPerKg;
this.initialLayerHeatContent = this.layerHeatContent;
this.layerHeatContent = heatPerKg * this.mass;
this.initialLayerHeatContentPerKg = this.layerHeatContentPerKg;
this.layerHeatAdded = 0;
averageTemperatureK = new Average( (int)( 1.0 / astro.ap.timeStepYears ) );
// this.setLayerPhase();
this.setTickYears( tickYears );
// System.out.println( "tickYears " + tickYears + " tick duration " +this.tickDuration + " tick count " + this.tickResetCount );
// System.out.println( "Layer constructor current temperature " + this.currentTemperatureK + " heat content " + this.layerHeatContent );
// System.out.println( "Layer constructor tick duration " + this.tickDuration / astro.secondsPerYear + " count " + this.tickResetCount );
// phase change data
double mt0 = (int)this.materialType; // solid (ice)
double mt1 = mt0 + 0.1; // liquid (water)
double mt2 = mt0 + 0.2; // gas (water vapour)
this.material0 = new Material( astro, mt0 );
this.material1 = new Material( astro, mt1 );
this.material2 = new Material( astro, mt2 );
// double spHt0 = this.layerSpecificHeat; // use current specific heat
this.spHt0 = material0.specificHeat;
this.spHt1 = material1.specificHeat;
this.spHt2 = material2.specificHeat;
pt = new PhaseTemperature( this );
// set temperatures using heat content of layer
this.currentTemperatureK = pt.temperature( layerHeatContentPerKg );
this.nextTemperatureK = this.currentTemperatureK;
this.initialTemperature = this.currentTemperatureK;
// if ( (int)mt == 5 ) this.materialType = 5.0; // fix for .hcm bug?
// if ( mt == 0 ) System.out.println( "this.currentTemperatureK " + this.currentTemperatureK );
}
// copy the characteristics of another layer
void copyLayer( ConductiveLayer layer ) {
// this.radius = layer.radius; not copied because it's a different layer
// nBand not copied because it's a different layer
// nBand is set by asteroid.reSynchronise();
this.materialType = layer.materialType;
this.material = layer.material;
this.length = layer.length;
this.width = layer.width;
this.height = layer.height;
this.area = layer.area;
this.volume = layer.volume;
this.mass = layer.mass;
this.currentTemperatureK = layer.currentTemperatureK;
this.nextTemperatureK = layer.nextTemperatureK;
this.layerDensity = layer.layerDensity;
this.layerConductivity = layer.layerConductivity;
this.layerSpecificHeat = layer.layerSpecificHeat;
this.layerLatentHeatOfFusion = layer.layerLatentHeatOfFusion;
this.layerEmissivity = layer.layerEmissivity;
this.layerAbsorptivity = layer.layerAbsorptivity;
this.layerThermalCapacity = layer.layerThermalCapacity;
this.heatProduced = layer.heatProduced;
this.layerHeatContent = layer.layerHeatContent;
this.albedo = layer.albedo;
}
void setActiveCalculation( boolean astate ) {
this.activeCalculation = astate;
}
// get a unique number for the layer
int getUniqueLayerNumber() {
astro.uniqueNumber++;
this.uniqueNumber = astro.uniqueNumber;
return( astro.uniqueNumber );
}
void restoreInitialTemperature() {
this.currentTemperatureK = this.initialTemperature;
this.nextTemperatureK = this.initialTemperature;
this.layerHeatContentPerKg = this.initialLayerHeatContentPerKg;
this.layerHeatContent = this.initialLayerHeatContent;
}
void fixTemperature() {
this.fixTemperature = true;
}
void unfixTemperature() {
this.fixTemperature = false;
}
// Layers are thin slices of a sphere, with this.width * this width cross-sectional area at top.
double layerVolume( double radiusTopOfLayer, double radiusBottomOfLayer, double sliceWidth ) {
double v = sphereVolume( radiusTopOfLayer ) - sphereVolume( radiusBottomOfLayer );
double fraction = sliceWidth * sliceWidth / sphereSurfaceArea( radiusTopOfLayer );
return( fraction * v );
}
// find layer volume from mass and density
double layerVolume( double mass, double density ) {
double v = mass / density;
return( v );
}
// layer is part of a column with surface area of one square metre on asteroid of radius referenceRadius
double radiusTopOfLayer( double layerVolume, double radiusBottomOfLayer, double referenceRadius ) {
double asteroidSurfaceArea = sphereSurfaceArea( referenceRadius );
double columnFractionOfSphere = 1.0 / asteroidSurfaceArea;
double rCubed = Math.pow( radiusBottomOfLayer, 3.0 ) + ( 3.0 * layerVolume ) / ( 4.0 * columnFractionOfSphere * Math.PI );
double r = Math.pow( rCubed, ( 1.0/3.0 ) );
return( r );
}
// set layer height, width, and top surface area
void setLayerDimensions( double radiusTopOfLayer, double radiusBottomOfLayer, double referenceRadius ) {
this.height = radiusTopOfLayer - radiusBottomOfLayer;
this.width = ( radiusTopOfLayer / referenceRadius ); // 1 m at asteroid surface
this.area = this.width * this.width; // 1 square m at surface
}
double layerMass() {
return( this.volume * this.material.density );
}
double layerThermalCapacity() {
return( this.volume * this.material.density * this.material.specificHeat );
}
double heatGenerated() {
return( this.heatProduced );
}
// fractional part of floating point number
double fPart( double number ) {
return( number - iPart( number ) );
}
double layerResistance() {
this.layerResistance = this.height / ( this.layerConductivity * this.area );
return( this.layerResistance );
}
double sphereVolume( double r ) {
double volume = 1.3333333333 * Math.PI * r * r * r;
return( volume );
}
double sphereSurfaceArea( double r ) {
double a = 4.0 * Math.PI * r * r;
return( a );
}
/*
// layer heat content including latent heat
double heatContent( double temperature ) {
double heatContent = 0;
Material material0;
Material material1;
Material material2;
double mt0 = (int)this.materialType; // solid (ice)
double mt1 = mt0 + 0.1; // liquid (water)
double mt2 = mt0 + 0.2; // gas (water vapour)
material0 = new Material( astro, mt0 );
material1 = new Material( astro, mt1 );
material2 = new Material( astro, mt2 );
// double spHt0 = this.layerSpecificHeat; // use current specific heat
double spHt0 = material0.specificHeat;
double spHt1 = material1.specificHeat;
double spHt2 = material2.specificHeat;
if ( this.nBand == 72 ) System.out.println( "hc " + heatContent + " phasechange " + this.phaseChange );
if ( this.phaseChange ) {
// if phase change in process layer will be at either phasetemperature1 (fusion) or phasetemperature2 (vapourisation)
if ( this.currentTemperatureK == material0.phaseTemperature1 ) {
heatContent = this.mass * spHt0 * material0.phaseTemperature1;
heatContent += this.mass * material0.latentHeatOfFusion * Math.abs( material0.phaseTemperature1 - this.pseudoTemperatureK );
// System.out.println( "heatContent2 " + heatContent + " " + this.mass + " " +material0.latentHeatOfFusion + " " + Math.abs( material0.phaseTemperature1 - this.pseudoTemperatureK ));
} else if ( this.currentTemperatureK == material0.phaseTemperature2 ) {
heatContent = this.mass * spHt0 * material0.phaseTemperature1;
heatContent += this.mass * material0.latentHeatOfFusion;
heatContent += this.mass * spHt1 * (material1.phaseTemperature2 - material1.phaseTemperature1);
heatContent += this.mass * Math.abs( material1.phaseTemperature2 - this.pseudoTemperatureK ) * material1.latentHeatOfEvaporation;
}
} else {
// if between phases
if ( temperature > material0.phaseTemperature2 ) {
heatContent = this.mass * spHt0 * material0.phaseTemperature1;
heatContent += this.mass * material0.latentHeatOfFusion;
heatContent += this.mass * spHt1 * ( material1.phaseTemperature2 - material1.phaseTemperature1 );
heatContent += this.mass * material1.latentHeatOfEvaporation;
heatContent += this.mass * spHt2 * ( temperature - material2.phaseTemperature2 );
} else if ( temperature > material0.phaseTemperature1 ) {
heatContent = this.mass * spHt0 * material0.phaseTemperature1;
heatContent += this.mass * material0.latentHeatOfFusion;
heatContent += this.mass * spHt1 * ( temperature - material1.phaseTemperature1 );
} else {
heatContent = this.mass * spHt0 * temperature;
}
}
if ( this.nBand == 72 ) System.out.println( "hc " + heatContent + " phasechange " + this.phaseChange );
return( heatContent );
}
*/
// early version prior to jan 2019
int setLayerPhase() {
if ( material.isSolid( this.currentTemperatureK ) ) {
this.phase = 0;
} else if ( material.isLiquid( this.currentTemperatureK ) ) {
this.phase = 1;
} else {
this.phase = 2;
}
return( this.phase );
}
// jan 2019 set new phase and material type if necessary
double phaseCheck( double t) {
int newphase = getPhase( t );
if ( newphase > -1 ) {
if ( newphase != this.phase) {
this.phase = newphase;
this.materialType = getMaterialType( this.materialType, newphase );
System.out.println( (float)astro.ap.year + " band " + this.nBand + " new phase " + newphase + " new material " + this.materialType );
} else {
// keep current materialType unchanged
}
} else {
// keep current materialType unchnanged during phase change
}
return( this.materialType );
}
// jan 2019
int getPhase( double t ) {
int lphase = -2;
if ( t == this.material.phaseTemperature1 || t == this.material.phaseTemperature2 ) {
lphase = -1; // phase change in process
} else {
if ( t < this.material.phaseTemperature1 ) {
lphase = 0; // solid
} else if ( t > this.material.phaseTemperature2 ) {
lphase = 2; // vapour
} else {
lphase = 1; // liquid
}
}
return( lphase );
}
// generic material type is the integer number of the material
// i.e. materialType 5.2 is generic material type 5.
int genericMaterialType() {
return( (int)this.materialType );
}
// jan 2019
double getMaterialType( double mt, int mphase ) {
double type = iPart( mt ) + (double)mphase / 10.0;
return( type );
}
// interger part of floating point number
int iPart( double number ) {
return ( (int)number );
}
int setTickYears( double years ) {
if ( years == 0 ) {
// default setting 0 is to use same timestep as ap
this.tickDuration = astro.ap.timeStepYears * astro.secondsPerYear;
this.tickResetCount = 1;
resetTickCounter();
} else if ( years > 0 ) {
if ( years > 0.1 * this.layerTimeConstant / astro.secondsPerYear ) {
years = 0.1 * this.layerTimeConstant / astro.secondsPerYear;
}
this.tickResetCount = (int)( years / astro.ap.timeStepYears );
if ( this.tickResetCount < 0 ) this.tickResetCount = 1;
this.tickDuration = years * astro.secondsPerYear;
}
return( this.tickResetCount );
}
int setupTick( int c ) {
if ( c == 0 ) {
// default setting 0 is to use same timestep as ap
this.tickDuration = astro.ap.timeStepYears * astro.secondsPerYear;
this.tickResetCount = 1;
resetTickCounter();
} else if ( c > 0 ) {
this.tickResetCount = (int)(this.layerTimeConstant / this.tickDuration);
this.tickResetCount = this.tickResetCount / c;
if ( this.tickResetCount < 0 ) this.tickResetCount = 1;
this.tickDuration = astro.ap.timeStepYears * astro.secondsPerYear * this.tickResetCount;
}
return( this.tickResetCount );
}
void resetTickCounter() {
this.tickCounter = this.tickResetCount;
}
}
class Material{
GeologicalColumn astro;
String name;
double index;
double albedo; // negative albedo is transparent, use albedo of layer below
double density;
double conductivity;
double specificHeat;
double emissivity;
double absorptivity;
double latentHeatOfFusion;
double latentHeatOfEvaporation;
double volumetricFractionGranite = 0;
double radioactiveHeatGenerated;
double phaseTemperature1; // melting point
double phaseTemperature2; // boiling point
double expansivity; // coefficient of linear thermal expansion
double modulusElasticity;
int colourIndex[] = new int[3]; // 0 solid, 1 liquid, 2 vapaor
double[][] iceCharacteristics = {
// copied from https://www.engineeringtoolbox.com/ice-thermal-properties-d_576.html
// degC density k sp ht (kJ/kg)
{ 0, 916.2, 2.22, 2.050, },
{ -5, 917.5, 2.25, 2.027, },
{ -10, 918.9, 2.30, 2.000, },
{ -15, 919.4, 2.34, 1.972, },
{ -20, 919.4, 2.39, 1.943, },
{ -25, 919.6, 2.45, 1.913, },
{ -30, 920.0, 2.50, 1.882, },
{ -35, 920.4, 2.57, 1.851, },
{ -40, 920.8, 2.63, 1.818, },
{ -50, 921.6, 2.76, 1.751, },
{ -60, 922.4, 2.90, 1.681, },
{ -70, 923.3, 3.05, 1.609, },
{ -80, 924.1, 3.19, 1.536, },
{ -90, 924.9, 3.34, 1.463, },
{-100, 925.7, 3.48, 1.389, },
{-150, 926.7, 3.52, 1.360, }, // extrapolated value (v approx)
{-200, 927.7, 3.56, 1.340, }, // extrapolated value (v approx)
{-250, 928.7, 3.60, 1.320, }, // extrapolated value (v approx)
};
Material() { }
Material( GeologicalColumn a, double i, String name, double d, double c, double s, double lhf, double lhv, double fgranite, double hgranite ) {
astro = a;
this.name = name;
this.index = i;
this.density = d;
this.conductivity = c;
this.specificHeat = s;
this.latentHeatOfFusion = lhf;
this.latentHeatOfEvaporation = lhv;
this.volumetricFractionGranite = fgranite;
this.radioactiveHeatGenerated = hgranite;
}
// select from existing materials
Material( GeologicalColumn a, double i ) {
astro = a;
setMaterialCharacteristics(i, 0);
}
// select from existing materials
Material( double i ) {
setMaterialCharacteristics(i, 0);
}
void setMaterialCharacteristics( double i, int calledBy ) {
int multiple = (int)( i * 100.01 );
switch ( multiple ) {
case 0:
case 10:
case 20:
// thin surface layer has zero thickness and therefore zero mass volume etc
this.name = "surface";
this.index = i;
this.albedo = -1; //-1 means albedo is of material in layer below
this.density = 7870;
this.conductivity = 79.5; // iron
this.specificHeat = 460.5; // iron
this.phaseTemperature1 = 100000000; // no phase
this.phaseTemperature2 = 100000000; // no phase change
this.emissivity = 0.93; // ???
this.absorptivity = 0.93;
this.expansivity = 0;
this.modulusElasticity = -0;
this.colourIndex[0] = 0; // black
this.colourIndex[1] = 0; //
this.colourIndex[2] = 0;
break;
case 100:
case 110:
case 120:
// basalt rock
this.name = "basalt asthenosphere";
this.index = i;
this.albedo = 0.3;
this.density = 3000.0; // kg / m^3 basalt
this.conductivity = 2.20; // W/ m deg K basalt
this.specificHeat = 840; // J/ kg degK basalt
this.phaseTemperature1 = 1300;
this.latentHeatOfFusion = 418400; // 100 calories / gm
this.phaseTemperature2 = 7001;
this.latentHeatOfEvaporation = 840; // set same as specific heat
this.volumetricFractionGranite = 0;
this.radioactiveHeatGenerated = 1E-9; // modern terrestrial granite
this.emissivity = 0.95; // https://www.optotherm.com/emiss-table.htm
this.absorptivity = 0.95;
this.expansivity = 8.0e-6; // 10^-6 same as granite
this.modulusElasticity = 50E9; // http://community.dur.ac.uk/~des0www4/cal/dams/geol/mod.htm
this.colourIndex[0] = 0; // brown
this.colourIndex[1] = 5; // red
this.colourIndex[2] = 7; // yellow
break;
case 200:
case 220:
// granite rock
this.name = "granite crust";
this.index = i;
// this.albedo = 0.4;
this.albedo = 0.4;
this.density = 2650.0; // kg / m^3
this.conductivity = 2.12; // dry high porosity W/ m deg K
// this.conductivity = 3.98; // wet low porosity W/ m deg K (may be as high as 5.00)
// this.conductivity = 6.36; // 3 x standard conductivity)
this.specificHeat = 790; // J/ kg degK engineering toolbox
this.phaseTemperature1 = 1523; // 1250 °C granite melts at about 1215° to 1260° C
this.latentHeatOfFusion = 418000; // conversion of 100 calories/gm
this.phaseTemperature2 = 7001; // large number so it won't happen
this.latentHeatOfEvaporation = 790; // set same as specific heat
this.volumetricFractionGranite = 1;
this.radioactiveHeatGenerated = 1E-9; // modern terrestrial granite
this.emissivity = 0.95; // https://www.optotherm.com/emiss-table.htm
this.absorptivity = 0.95;
this.expansivity = 8.0e-6; // 10^-6
this.modulusElasticity = 40E9; //http://community.dur.ac.uk/~des0www4/cal/dams/geol/mod.htm
this.colourIndex[0] = 9; // brown
this.colourIndex[1] = 5; // red
this.colourIndex[2] = 7; // yellow
break;
case 205:
// foamed granite rock (aka pumice) copied from material index 11
// pumice (all the characteristics of granit except density and conductivity)
this.name = "pumice, or bubbly granite";
this.index = i;
this.albedo = 0.4;
this.density = 50.0; // slightly higher than polystyrene foam kg / m^3
this.conductivity = 0.02; // conductivity of polystyrene foam W/ m deg K
this.specificHeat = 790; // J/ kg degK engineering toolbox
this.phaseTemperature1 = 1523;
this.latentHeatOfFusion = 418000; // conversion of 100 calories/gm
this.phaseTemperature2 = 7001; // large number so it won't happen
this.latentHeatOfEvaporation = 460; // set same as specific heat
this.volumetricFractionGranite = 1;
this.radioactiveHeatGenerated = 1E-9; // modern terrestrial granite
this.emissivity = 0.95; // https://www.optotherm.com/emiss-table.htm
this.absorptivity = 0.95;
this.expansivity = 8.0e-6; // 10^-6
this.modulusElasticity = 40E9; //http://community.dur.ac.uk/~des0www4/cal/dams/geol/mod.htm
this.colourIndex[0] = 6; // brown
this.colourIndex[1] = 6; // red
this.colourIndex[2] = 7; // yellow
break;
case 210:
// molten granite rock is supposed to have the same characteristic as solid rock
this.name = "granite lava";
this.index = i;
this.albedo = 0.4;
// this.albedo = 0.50;
this.density = 2650.0; // kg / m^3
this.conductivity = 2.12; // dry high porosity W/ m deg K
// this.conductivity = 3.98; // wet low porosity W/ m deg K
this.specificHeat = 790; // J/ kg degK engineering toolbox
this.phaseTemperature1 = 1523; // 1250 °C granite melts at about 1215° to 1260° C
this.latentHeatOfFusion = 418000; // conversion of 100 calories/gm
this.phaseTemperature2 = 7001; // large number so it won't happen
this.latentHeatOfEvaporation = 460; // set same as specific heat
this.volumetricFractionGranite = 1;
this.radioactiveHeatGenerated = 1E-9; // modern terrestrial granite
this.emissivity = 0.95; // https://www.optotherm.com/emiss-table.htm
this.absorptivity = 0.95;
this.expansivity = 8.0e-6; // 10^-6
this.modulusElasticity = 40E9; //http://community.dur.ac.uk/~des0www4/cal/dams/geol/mod.htm
this.colourIndex[0] = 9; // brown
this.colourIndex[1] = 5; // red
this.colourIndex[2] = 7; // yellow
break;
case 300:
// water ice
this.name = "ice";
this.index = 3.0;
this.albedo = 0.6; // less than snow
this.density = 917; // density of ice kg/ m^3
this.conductivity = 2.18; // conductivity of Ice W / m degK
this.specificHeat = 1389; // specific heat of ice, J / kg deg K
this.phaseTemperature1 = 273;
this.latentHeatOfFusion = 334000; // J/kg
this.phaseTemperature2 = 373;
this.latentHeatOfEvaporation = 2264705; // J/kg
this.volumetricFractionGranite = i - (int)i;
this.radioactiveHeatGenerated = 0;
this.emissivity = 0.85; // snow 0.8 - 0.9
this.absorptivity = 0.85;
this.expansivity = 51.0;
this.modulusElasticity = 0; // unknown
this.colourIndex[0] = 3; // mid gray (ice)
this.colourIndex[1] = 12; // cyan (water)
this.colourIndex[2] = 7; // yellow (steam)
break;
case 305:
// snow, which is a form of ice, but less dense and less conductive
this.name = "snow";
this.index = 3.05;
// this.albedo = 0.8; // fresh snow albedo
// this.albedo = 0.80; // fresh snow albedo
this.albedo = 0.80;
this.density = 30; // density of dry light snow kg/ m^3
this.conductivity = 0.02; // conductivity of snow W / m degK ( same as exp polystyrene)
// this.conductivity = 0.006; // conductivity of PP snow 103 kg/m3 https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2011GL049234
// this.conductivity = 0.08; // conductivity of fresh snow W / m degK
this.specificHeat = 2025; // specific heat of ice, J / kg deg K
this.phaseTemperature1 = 273;
this.latentHeatOfFusion = 334000; // J/kg
this.phaseTemperature2 = 373;
this.latentHeatOfEvaporation = 2264705; // J/kg
this.volumetricFractionGranite = i - (int)i;
this.radioactiveHeatGenerated = 0;
this.emissivity = 0.85; // snow 0.8 - 0.9
this.absorptivity = 0.85;
this.expansivity = 51.0;
this.modulusElasticity = 0; // unknown
this.colourIndex[0] = 4; // mid gray (ice)
this.colourIndex[1] = 12; // cyan (water)
this.colourIndex[2] = 7; // yellow (steam)
break;
case 310:
// liquid water
this.name = "water";
this.index = 3.1;
this.albedo = -0.3; // negative means transparent
this.density = 999.8; // water density at 273 K
// this.conductivity = 0.59; // water conductivity at 273 K
this.conductivity = 0.56; // temporary water conductivity to emulate convection
this.specificHeat = 4217; // water specific heat at 273 K
this.phaseTemperature1 = 273;
this.latentHeatOfFusion = 334000; // J/kg
this.phaseTemperature2 = 373;
this.latentHeatOfEvaporation = 2264705; // J/kg
this.volumetricFractionGranite = i - (int)i;
this.radioactiveHeatGenerated = 0;
this.emissivity = 0.96; // water, snow 0.8 - 0.9
this.absorptivity = 0.96;
this.expansivity = 0; // water dees not expand much
this.modulusElasticity = 2.2E9; // 2.2 GN/m2,
this.colourIndex[0] = 4; // mid gray
this.colourIndex[1] = 12; // blue
this.colourIndex[2] = 7; // light gray
break;
case 320:
// steam or water vapour
// using characteristics f water for now
this.name = "steam";
this.index = 3.2;
this.albedo = -0.3; // transparent
this.density = 999.8; // water density at 273 K
this.conductivity = 0.59; // water conductivity at 273 K
this.specificHeat = 1996; // water specific heat at 273 K
this.phaseTemperature1 = 273;
this.latentHeatOfFusion = 334000; // J/kg
this.phaseTemperature2 = 373;
this.latentHeatOfEvaporation = 2264705; // J/kg
this.volumetricFractionGranite = i - (int)i;
this.radioactiveHeatGenerated = 0;
this.emissivity = 0.96; // ???
this.absorptivity = 0.96;
this.expansivity = 0; // no idea!
this.modulusElasticity = -0; // unknown
this.colourIndex[0] = 4; // mid gray
this.colourIndex[1] = 12; // cyan
this.colourIndex[2] = 1; // light gray
break;
case 500:
case 510:
case 520:
// air
this.name = "air";
this.index = i;
this.albedo = -0.3;
this.density = 1.293; // density of air kg/ m^3
// this.density = 2.586; // double density of air kg/ m^3
this.conductivity = 0.0243; // conductivity of air W / m degK
this.specificHeat = 1000; // specific heat of air, J / kg deg K
this.phaseTemperature1 = 58;
this.latentHeatOfFusion = 1000; // set same as specific heat
this.phaseTemperature2 = 76.3;
this.latentHeatOfEvaporation = 1000; // set same as specific heat
this.volumetricFractionGranite = i - (int)i;
this.radioactiveHeatGenerated = 0;
// this.emissivity = 0.958; // RESULTS IN 300 k AIR TEMPERATURES AT EQUATOR
// this.absorptivity = 0.958;
// this.emissivity = 0.95;
// this.absorptivity = 0.95;
// this.emissivity = 0.999; // effect of greenhouse gases?
// this.absorptivity = 0.999;
// this.emissivity = 0.80; // ACS value
// this.absorptivity = 0.80;
this.emissivity = 0.91; // prelininaries 0.90 -> 288.08522217136186 K using Solar=1, 0.91 -> 288.04 using Solar=2.
this.absorptivity = 0.91;
// this.emissivity = 0.911;
// this.absorptivity = 0.911;
this.expansivity = 0.00369; // 1/K coefficient at 0°C and 1 bar https://www.engineeringtoolbox.com/air-properties-d_156.html
this.modulusElasticity = 1.01325E5; // https://www.engineeringtoolbox.com/air-properties-d_156.html
this.colourIndex[0] = 10; // cyan
this.colourIndex[1] = 10; // cyan
this.colourIndex[2] = 10; // cyan ay
break;
case 600:
case 610:
case 620:
// planet inner core (iron?)
this.name = "inner iron core";
this.index = i;
this.albedo = 0.3;
this.density = 13000; // http://hyperphysics.phy-astr.gsu.edu/hbase/Geophys/earthstruct.html
this.conductivity = 79.5; // iron
this.specificHeat = 800.0; // file:///C:/Users/Me/Downloads/m00024-dat-914-core-properties.pdf
this.phaseTemperature1 = 6501; // ?
this.latentHeatOfFusion = 460.5; // set same as specific heat
this.phaseTemperature2 = 7001; // value to ensure doesn't vaporize
this.latentHeatOfEvaporation = 460.5; // set same as specific heat
this.volumetricFractionGranite = i - (int)i;
this.radioactiveHeatGenerated = 0;
this.emissivity = 0.87;
this.absorptivity = 0.87;
this.expansivity = 10.5E-6;
this.modulusElasticity = -0;
this.colourIndex[0] = 0; // black
this.colourIndex[1] = 6; // orange
this.colourIndex[2] = 7; // yellow
break;
case 700:
case 710:
case 720:
// planet outer core (iron?)
this.name = "outer iron core";
this.index = i;
this.albedo = 0.3;
this.density = 12200; // http://hyperphysics.phy-astr.gsu.edu/hbase/Geophys/earthstruct.html
this.conductivity = 79.5; // iron
this.specificHeat = 800.0; // file:///C:/Users/Me/Downloads/m00024-dat-914-core-properties.pdf
this.phaseTemperature1 = 6001; // ?
this.latentHeatOfFusion = 460.5; // set same as specific heat
this.phaseTemperature2 = 7001; // value to ensure doesn't vaporize
this.latentHeatOfEvaporation = 460.5; // set same as specific heat
this.volumetricFractionGranite = i - (int)i;
this.radioactiveHeatGenerated = 0;
this.emissivity = 0.87;
this.absorptivity = 0.87;
this.expansivity = 10.5E-6; // same as inner core
this.modulusElasticity = -0; // unknown
this.colourIndex[0] = 0; // black
this.colourIndex[1] = 6; // orange
this.colourIndex[2] = 7; // yellow
break;
case 800:
case 810:
case 820:
// planet lower mantle)
this.name = "lower mantle";
this.index = i;
this.albedo = 0.3;
this.density = 4600; // http://hyperphysics.phy-astr.gsu.edu/hbase/Geophys/earthstruct.html
this.conductivity = 6.0; // https://escholarship.org/uc/item/6ns1k73c
this.specificHeat = 1000.0; // https://www.colorado.edu/geolsci/courses/GEOL5700-PCE/convection.pdf
this.phaseTemperature1 = 5001; // ?
this.latentHeatOfFusion = 460.5; // set same as specific heat
this.phaseTemperature2 = 7001; // value to ensure doesn't vaporize
this.latentHeatOfEvaporation = 460.5; // set same as specific heat
this.volumetricFractionGranite = i - (int)i;
this.radioactiveHeatGenerated = 0;
this.emissivity = 0.87;
this.absorptivity = 0.87;
this.expansivity = 15.0E-6; // MgO
this.modulusElasticity = -0; // unknown
this.colourIndex[0] = 0; // black
this.colourIndex[1] = 6; // orange
this.colourIndex[2] = 7; // yellow
break;
case 900:
case 910:
case 920:
// planet upper mantle)
this.name = "upper mantle";
this.index = i;
this.albedo = 0.3;
this.density = 3600; // http://hyperphysics.phy-astr.gsu.edu/hbase/Geophys/earthstruct.html
this.conductivity = 4.5; // guessed
this.specificHeat = 1000.0; // https://www.colorado.edu/geolsci/courses/GEOL5700-PCE/convection.pdf
this.phaseTemperature1 = 5001; // ?
this.latentHeatOfFusion = 460.5; // set same as specific heat
this.phaseTemperature2 = 7001; // value to ensure doesn't vaporize
this.latentHeatOfEvaporation = 460.5; // set same as specific heat
this.volumetricFractionGranite = i - (int)i;
this.radioactiveHeatGenerated = 0;
this.emissivity = 0.87;
this.absorptivity = 0.87;
this.expansivity = 30.0E-6; // MgO
this.modulusElasticity = -0; // unknown
this.colourIndex[0] = 0; // black
this.colourIndex[1] = 6; // orange
this.colourIndex[2] = 7; // yellow
break;
case 1000:
case 1010:
case 1020:
// outer space (using characteristics of air)
this.name = "space";
this.index = i;
this.albedo = -0.3;
this.density = 0.00001; // density of air kg/ m^3
this.conductivity = 0.00001; // conductivity of air W / m degK
this.specificHeat = 0.0001; // specific heat of air, J / kg deg K
this.phaseTemperature1 = 58;
this.latentHeatOfFusion = 1000; // set same as specific heat
this.phaseTemperature2 = 76.3;
this.latentHeatOfEvaporation = 1000; // set same as specific heat
this.volumetricFractionGranite = i - (int)i;
this.radioactiveHeatGenerated = 0;
this.emissivity = 0.95; // ???
this.absorptivity = 0.95;
this.expansivity = 0.00369; // 1/K coefficient at 0°C and 1 bar https://www.engineeringtoolbox.com/air-properties-d_156.html
this.modulusElasticity = 1.01325E5; // https://www.engineeringtoolbox.com/air-properties-d_156.html
this.colourIndex[0] = 10; // mid gray
this.colourIndex[1] = 10; // cyan
this.colourIndex[2] = 10; // light gray
break;
case 1100:
case 1110:
case 1120:
// pumice (all the characteristics of granit except density and conductivity)
this.name = "pumice, or bubbly granite";
this.index = i;
this.albedo = 0.3;
this.density = 50.0; // slightly higher than polystyrene foam kg / m^3
this.conductivity = 0.02; // conductivity of polystyrene foam W/ m deg K
this.specificHeat = 790; // J/ kg degK engineering toolbox
this.phaseTemperature1 = 1250;
this.latentHeatOfFusion = 418000; // conversion of 100 calories/gm
this.phaseTemperature2 = 10000; // large number so it won't happen
this.latentHeatOfEvaporation = 460; // set same as specific heat
this.volumetricFractionGranite = 1;
this.radioactiveHeatGenerated = 1E-9; // modern terrestrial granite
this.emissivity = 0.95; // https://www.optotherm.com/emiss-table.htm
this.absorptivity = 0.95;
this.expansivity = 8.0e-6; // 10^-6
this.modulusElasticity = 40E9; //http://community.dur.ac.uk/~des0www4/cal/dams/geol/mod.htm
this.colourIndex[0] = 6; // brown
this.colourIndex[1] = 6; // red
this.colourIndex[2] = 7; // yellow
break;
case 1200:
case 1210:
case 1220:
// granite rock
this.name = "chalk";
this.index = i;
this.albedo = 0.7;
this.density = 1790.0; // kg / m^3 nora.nerc.ac.uk/7752/
// this.conductivity = 1.78; // W/ m deg K https://qjegh.lyellcollection.org/content/early/2018/05/28/qjegh2017-127
this.conductivity = 0.09; // wet low porosity W/ m deg K (may be as high as 5.00)
this.specificHeat = 750; // J/ kg degK https://www.engineeringtoolbox.com/specific-heat-capacity-d_391.html
this.phaseTemperature1 = 1612; // https://en.wikipedia.org/wiki/Calcium_carbonate
this.phaseTemperature1 = 7000; // guessed
this.latentHeatOfFusion = 418000; // guessed same as granite
this.phaseTemperature2 = 7001; // guessed large number so it won't happen
this.latentHeatOfEvaporation = 720; // guessed set same as specific heat
this.volumetricFractionGranite = 0;
this.radioactiveHeatGenerated = 0;
this.emissivity = 0.95; // https://www.optotherm.com/emiss-table.htm
this.absorptivity = 0.95;
this.expansivity = 8.0e-6; // 10^-6 1/K same as granite
this.modulusElasticity = 40E9; // same as granite
this.colourIndex[0] = 4; // white
this.colourIndex[1] = 12; // cyan (water)
this.colourIndex[2] = 7; // yellow (steam)
break;
default:
System.out.println( "Can't find material " + i );
break;
}
}
//
// for use when ice (index 3.0) is combined with water (index 3.2)
boolean isWater( double t ) {
boolean isWater = false;
if ( this.index == 3 ) {
// water phase change over 1 degree K
if ( t > this.phaseTemperature1+0.01 && t < this.phaseTemperature2 ) {
isWater = true;
}
}
return ( isWater );
}
// for use when ice (material index 3) is combined with water (material index 4)
boolean isIce( double t ) {
boolean isIce = false;
if ( this.index == 3 ) {
// water phase change over 1 degree K
if ( t < this.phaseTemperature1+0.01 && t < this.phaseTemperature2 ) {
isIce = true;
}
}
return ( isIce );
}
boolean isLiquid( double t ) {
boolean isLiquid = false;
if ( t > this.phaseTemperature1+0.01 && t < this.phaseTemperature2 ) {
isLiquid = true;
}
return ( isLiquid );
}
boolean isSolid( double t ) {
boolean isSolid = false;
if ( t < this.phaseTemperature1+0.01 && t < this.phaseTemperature2 ) {
isSolid = true;
}
return ( isSolid );
}
boolean isAir( double t ) {
boolean isAir = false;
if ( this.index == 5 ) {
isAir = true;
}
return ( isAir );
}
int indexIceTemperature( double degK ) {
int index = 0;
for ( int n=0; n<18; n++ ) {
if ( celsiusToKelvin( iceCharacteristics[n][0] ) < degK ) {
index = n;
break;
}
}
return( index);
}
double celsiusToKelvin( double degC ) {
return( degC + 273 );
}
}
class PhaseTemperature{
ConductiveLayer layer;
double t;
double phase;
double spHt0; // solid specific heat
double latHt1; // latenet heat of fusion
double spHt1; // liquid specific heat
double latHt2; // latenet heat of evaporation
double spHt2; // vapour specific heat
double phaseT1; // solid-liquid phase change temperature
double phaseT2; // liquid-vapour phase change temperature
double phaseT3; // non-existent high temperature
double phaseStep[][] = new double[6][5];
double FractionPhaseChange;
// default constructor
PhaseTemperature() {
}
// find heat content at each phase step
PhaseTemperature( ConductiveLayer layer ) {
this.layer = layer;
double mt0 = layer.materialType;
double mt1 = (double)(int)mt0 + 0.1;
Material material1 = new Material( mt1 );
double mt2 = (double)(int)mt0 + 0.2;
Material material2 = new Material( mt2 );
this.spHt0 = this.layer.material.specificHeat; // layer solid phase specific heat
this.spHt1 = material1.specificHeat; // liquid phase specific heat
this.spHt2 = material2.specificHeat; // vapour phase sepcific heat
this.latHt1 = this.layer.material.latentHeatOfFusion;
this.latHt2 = this.layer.material.latentHeatOfEvaporation;
this.phaseT1 = this.layer.material.phaseTemperature1;
this.phaseT2 = this.layer.material.phaseTemperature2;
// System.out.println( "material type " + mt0 + " spht0 " + this.spHt0 );
// System.out.println( "material type " + mt1 + " spht1 " + this.spHt1 );
// System.out.println( "material type " + mt2 + " spht2 " + this.spHt2 );
this.phaseT3 = this.phaseT2 + this.spHt2; // imaginary temperature
// absolute zero
phaseStep[0][0] = 0; // temperature change K
phaseStep[0][1] = 0; // temperature K
phaseStep[0][2] = 0; // specific or latent heat
phaseStep[0][3] = 0; // heat content
// solid phase
phaseStep[1][0] = this.phaseT1; // temperature change K
phaseStep[1][1] = this.phaseT1; // temperature K
phaseStep[1][2] = this.spHt0 * this.phaseT1; // specific or latent heat
phaseStep[1][3] = phaseStep[1][2]; // heat content
// phase change solid to liquid
phaseStep[2][0] = 0; // temperature change K
phaseStep[2][1] = this.phaseT1; // temperature K
phaseStep[2][2] = this.latHt1; // specific or latent heat
phaseStep[2][3] = phaseStep[1][3] + this.latHt1; // heat content
// liquid phase
phaseStep[3][0] = this.phaseT2 - this.phaseT1; // temperature change K
phaseStep[3][1] = this.phaseT2; // temperature K
phaseStep[3][2] = phaseStep[3][0] * this.spHt1; // specific or latent heat
phaseStep[3][3] = phaseStep[2][3] + phaseStep[3][2]; // heat content
// phase change liquid to vapour
phaseStep[4][0] = 0; // temperature change K
phaseStep[4][1] = this.phaseT2; // temperature K
phaseStep[4][2] = this.latHt2; // specific or latent heat
phaseStep[4][3] = phaseStep[3][3] + this.latHt2; // heat content
// vapour phase
phaseStep[5][0] = this.phaseT3 - this.phaseT2; // temperature change K
phaseStep[5][1] = this.phaseT3; // temperature K
phaseStep[5][2] = phaseStep[5][0] * this.spHt1; // specific or latent heat
phaseStep[5][3] = phaseStep[4][3] + phaseStep[5][2]; // heat content
}
// find temperature and phase from layer heat content q
double temperature( double layerHeatContentPerKg ) {
int n =0;
this.t = 0;
this.phase = 0;
double q = layerHeatContentPerKg;
for ( n=0; n<=5; n++ ) {
if ( q < phaseStep[n][3] ) {
this.phase = (double)(n - 1) / 2.0;
// where t = temperature K and q = heat J
// t = t[n] + ( q - q[n-1] ) * T step change / q step change
this.t = phaseStep[n-1][1] + ( q - phaseStep[n-1][3] ) * phaseStep[n][0] / phaseStep[n][2] ;
this.FractionPhaseChange = ( q - phaseStep[n-1][3] ) / phaseStep[n][2] ;
break;
}
}
layer.nextTemperatureK = this.t;
if ( layer.materialType == 3.05 && phase != layer.nextPhase ) {
// System.out.println( layer.astro.ap.year + " snow phase change to " +this.phase );
layer.astro.ap.snowHeatAdded = 0;
}
layer.nextPhase = this.phase;
return( this.t );
}
// find layer heat content from temperature (added sep 2019)
double heatContentPerKg( double layerTemperature ) {
double q = 0;
for ( int n=0; n<=5; n++ ) {
if ( layerTemperature < phaseStep[n][1] ) {
q = phaseStep[n][3] + ( layerTemperature - phaseStep[n][1] ) * phaseStep[n][2] / phaseStep[n][0];
break;
}
}
return( q );
}
}
/*********** moving graph canvas *****************************************/
class ItemGraphCanvas extends Canvas {
Glaciation ap;
Dimension fielddimension;
// Image fieldImage;
// Graphics field2Graphics;
int day;
int maxstreams = 6;
int nstreams;
double lasty[] = new double[maxstreams];
double nexty[] = new double[maxstreams];
String streamLabel[] = new String[maxstreams];
double streamMaxValue[] = new double[maxstreams];
double streamMinValue[] = new double[maxstreams];
boolean plot[] = new boolean[maxstreams];
boolean overlappingStreams = true;
// int graph_height = 100;
int graphymax = 100;
boolean started;
Font font;
FontMetrics fm;
int redLetterDay;
String ylabel = new String();
Dimension d;
int sideband = 25;
int lowband = 15;
long calendar = 0; // count keeping elapsed time
double calendarStep; // actual duration of each pixel calendar step (e.g. years)
long calendarDisplayInterval;
String calendarLabel; // label to attach to calanddar dates
boolean autoCalendar = true; // if true, created dates automatically
Color colorpalette[] = new Color[15];
Image historicalRecord; // rolling image of historical graphs
Graphics historicalGraphics;
Dimension hr; // width and height of historical record
Image currentWindow; // current window onto historical record
Graphics currentGraphics;
Dimension cw; // width and height of historical record
// Constructor allows access to Field objects by this class
public ItemGraphCanvas( Glaciation app, Graphics g, int rollingGraphWidth, int rollingGraphHeight, int ns, boolean overlap ) {
super();
ap = app;
d = new Dimension( rollingGraphWidth, rollingGraphHeight );
// graph_height = d.height;
nstreams = ns;
overlappingStreams = overlap;
for ( int n = 0; n < nstreams; n++ ) {
lasty[n] = -100.0; nexty[n] = -100.0; streamLabel[n] = " ";
plot[n] = true;
}
ylabel = "500";
day = 0;
started = false;
redLetterDay = -1;
font = new Font( "Helvetica", Font.PLAIN, 12 );
calendarDisplayInterval = 100000;
calendarLabel = "e5";
colorpalette[0] = new Color(000, 000, 000); // black
colorpalette[1] = new Color( 64, 64, 64); // light grey
colorpalette[2] = new Color(128, 128, 128); // mid gray
colorpalette[3] = new Color(192, 192, 192); // darke gray
colorpalette[4] = new Color(255, 255, 255); // white
colorpalette[5] = Color.RED;
colorpalette[6] = Color.ORANGE;
colorpalette[7] = Color.YELLOW;
colorpalette[8] = Color.CYAN;
colorpalette[9] = new Color( 85, 52, 52); // brown?
colorpalette[10] = new Color(135,206, 250); // sky blue
colorpalette[11] = new Color(240,240, 240); // steam
colorpalette[12] = new Color( 0, 0, 255); // blue
colorpalette[13] = new Color( 0, 0, 255); // blue
colorpalette[14] = new Color( 0, 0, 255); // blue
}
void setNumberOfStreams( int n ) {
nstreams = n;
}
public void addStream() {
nstreams++;
}
public void removeStream() {
nstreams--;
}
public void setStreamLabel( int stream, String s ) {
streamLabel[stream] = s;
}
// set max and min values that are to be visible
public void setStreamMaxMin( int stream, double vmax, double vmin ) {
streamMaxValue[ stream ] = vmax;
streamMinValue[ stream ] = vmin;
}
public void setYscale( int ymax ) {
graphymax = ymax;
ylabel = Integer.toString( ymax );
}
public void stop() {
// field2Graphics = null;
// fieldImage = null;
}
// public void redraw() { repaint(); }
// update calls paint(), otherwise clears the screen.
public void update(Graphics g) {
paint(g);
}
public void paint(Graphics fieldGraphics ) {
int n, h;
String s;
calendar += calendarStep;
// System.out.println( "year " + (float)ap.year + " calendar " + (float)calendar + " calendarDisplayInterval " + (float)calendarDisplayInterval);
// Create the offscreen graphics context, if no good one exists.
if ( (fieldGraphics == null) ) {
// fieldImage = createImage(d.width, d.height + lowband);
// fieldGraphics = fieldImage.getGraphics();
fieldGraphics.setColor(Color.WHITE);
fieldGraphics.fillRect(0, 0, d.width, d.height + lowband);
fieldGraphics.setFont( font );
// g.drawImage(fieldImage, 0, 0, this);
}
fm = fieldGraphics.getFontMetrics( font );
if ( started ) {
// move image 1 pixel leftward
fieldGraphics.setPaintMode();
fieldGraphics.copyArea(1,0, d.width-sideband, d.height+lowband-1, -1, 0 );
// rub out unmoved vertical line
fieldGraphics.setColor(Color.WHITE);
fieldGraphics.fillRect(d.width-sideband-1, 0, d.width-1, d.height+lowband);
// add righthand text
fieldGraphics.setColor(Color.darkGray);
// fieldGraphics.drawString( ylabel, d.width-sideband+1, 10 );
// fieldGraphics.drawString( "0", d.width-sideband/3-1, d.height-lowband );
// draw horizontal gridlines
if ( this.overlappingStreams ) {
h = 0;
for ( n=0; n<=nstreams; n++ ) {
h = h + (d.height - lowband)/5;
fieldGraphics.setColor(Color.darkGray);
fieldGraphics.drawLine( d.width-sideband-2, h, d.width-sideband-1, h );
s = Integer.toString( (int)( streamMaxValue[0] ) );
fieldGraphics.drawString( s, d.width-sideband, rescale( 0, streamMaxValue[0] ) +9 );
s = Integer.toString( (int)( streamMinValue[0] ) );
fieldGraphics.drawString( s, d.width-sideband, rescale( 0, streamMinValue[0] ) -2 );
}
} else {
for ( n=0; n<nstreams; n++ ) {
fieldGraphics.setColor(Color.darkGray);
fieldGraphics.drawLine( d.width-sideband-2, rescale(n, 0), d.width-sideband-1, rescale(n, 0) );
s = Integer.toString( (int)( streamMaxValue[n] ) );
fieldGraphics.drawString( s, d.width-sideband, rescale( n, streamMaxValue[n] ) +9 );
s = Integer.toString( (int)( streamMinValue[n] ) );
fieldGraphics.drawString( s, d.width-sideband, rescale( n, streamMinValue[n] ) -2 );
}
}
// draw vertical gridlines and day numbers
if ( autoCalendar ) {
if ( calendar % calendarDisplayInterval == 0) {
// System.out.println( calendar + " calendar " + calendarDisplayInterval );
fieldGraphics.setColor(Color.darkGray);
fieldGraphics.drawLine( d.width-sideband-1, 0, d.width-sideband-1, d.height-lowband );
s = Integer.toString( (int)( calendar / calendarDisplayInterval ) ) + calendarLabel;
fieldGraphics.drawString( s, d.width-sideband-fm.stringWidth(s)-1, d.height-2 );
}
} else if ( (day % calendarDisplayInterval ) == 0 ) {
fieldGraphics.setColor(Color.darkGray);
fieldGraphics.drawLine( d.width-sideband-1, 0, d.width-sideband-1, d.height );
s = Integer.toString(day);
fieldGraphics.setColor(Color.darkGray);
fieldGraphics.drawString( s, d.width-sideband-fm.stringWidth(s)-1, d.height+lowband-1 );
}
// draw red vertical line to show event
if ( redLetterDay >= 0 ) {
fieldGraphics.setColor(Color.red);
fieldGraphics.drawLine( d.width-sideband-1, 0, d.width-sideband-1, d.height );
s = Integer.toString(redLetterDay);
fieldGraphics.drawString( s, d.width-sideband-fm.stringWidth(s)-3, fm.getAscent() + 1 );
}
// System.out.println( " overlap " + nstreams );
// draw graph lines and labels at right edge of graph
for ( n = nstreams-1; n >= 0; n-- ) {
if ( plot[n] ) {
fieldGraphics.setColor(Color.orange);
if ( n==0 ) fieldGraphics.setColor(Color.BLACK);
if ( n==1 ) fieldGraphics.setColor(Color.red);
if ( n==2 ) fieldGraphics.setColor(Color.green);
if ( n==3 ) fieldGraphics.setColor(Color.cyan);
if ( n==4 ) fieldGraphics.setColor(Color.orange);
if ( n==5 ) fieldGraphics.setColor(Color.darkGray);
fieldGraphics.drawLine( d.width-sideband-2, rescale( n, lasty[n] ), d.width-sideband-1, rescale( n, nexty[n] ) );
fieldGraphics.drawString( streamLabel[n], d.width-sideband+2, rescale( n, nexty[n] ) );
}
}
// fieldGraphics.setColor(Color.red);
// fieldGraphics.fillRect( 20, 20, 20, 20);
//Paint the image onto the screen.
// g.drawImage(fieldImage, 0, 0, this);
} else {
started = true;
fieldGraphics.setColor(Color.WHITE);
fieldGraphics.fillRect(0, 0, d.width, d.height);
}
for ( n = 0; n < nstreams; n++ ) lasty[n] = nexty[n];
}
// rescale y value
int rescale( double yvalue ) {
int rsv;
double gph = (double)d.height;
double yv = (double)yvalue;
double gymax = (double)graphymax;
double yl = gph - yv * ( gph / gymax );
rsv = (int)yl;
return( rsv );
}
// new 2018 rescale y value
// new code can handle +/- yvalue range
int rescale( int s, double yvalue ) {
int rsv;
double height = d.height - lowband;
if ( this.overlappingStreams ) {
double vscale = height / (streamMaxValue[s] - streamMinValue[s]);
double xaxis = vscale * streamMaxValue[s];
rsv = (int)( xaxis - yvalue * vscale );
// if ( s==0 ) System.out.println( "yvalue " + (float)yvalue + " vscale " + (float)vscale + " xaxis " + (float)xaxis + " rsv " + rsv);
} else {
// each stream height = ( d.height / (double)nstreams )
double vscale = height / ( (streamMaxValue[s] - streamMinValue[s]) * (double)nstreams );
double xaxis = ( (double)s * height / (double)nstreams ) + ( vscale * streamMaxValue[s] );
rsv = (int)( xaxis - yvalue * vscale );
}
return( rsv );
}
}
class LayerDisplay {
Glaciation ap;
int xcentre;
int ycentre;
int xtopleft;
int ytopleft;
int hr;
double jscale;
int xmax;
int ymax;
Color colorpalette[] = new Color[15];
boolean showLayerBoundaries = true;
double topLayerDepth = 40000;
int NBAND = 0;
int MATERIALTYPE = 1;
int RADIUS = 2;
int HEIGHT = 3;
int TEMPERATURE = 4;
int PHASE = 5;
int UHFLOWRATE = 6;
public LayerDisplay( Glaciation a ) {
this.ap = a;
colorpalette[0] = new Color(000, 000, 000); // black
colorpalette[1] = new Color( 64, 64, 64); // light grey
colorpalette[2] = new Color(128, 128, 128); // mid gray
colorpalette[3] = new Color(192, 192, 192); // darke gray
colorpalette[4] = new Color(255, 255, 255); // white
colorpalette[5] = Color.RED;
colorpalette[6] = Color.ORANGE;
colorpalette[7] = Color.YELLOW;
colorpalette[8] = Color.CYAN;
colorpalette[9] = new Color( 85, 52, 52); // brown?
colorpalette[10] = new Color(135,206, 250); // sky blue
colorpalette[11] = new Color(240,240, 240); // steam
colorpalette[12] = new Color( 0, 0, 255); // blue
colorpalette[13] = new Color( 0, 0, 255); // blue
colorpalette[14] = new Color( 0, 0, 255); // blue
}
void setDisplaySize( int x, int y ) {
xmax = x;
ymax = y;
// System.out.println( "layer xmax " + xmax + " ymax " + ymax );
}
void paint( Graphics g, GeologicalColumn asteroid, double topRockRadius ) {
double r;
// show 200 km of top layers of body
// double topLayerDepth = 200000.0;
double jscale = xmax / topLayerDepth; // screen to show top 40 km
jscale = 20.0 * jscale;
double tscale = xmax / 5200.0; // temperature scale
// clear window
g.setColor( Color.white );
g.fillRect( 0, 0, xmax, ymax );
// circle centre (pixels)
xcentre = xmax / 2;
ycentre = (int)( topRockRadius * jscale ) + ( ymax / 2 );
// paint atmosphere filled circle
// surface rock radius
int sr = (int)( ap.column[ ap.nviewLatitude ].graphicLayer[0][RADIUS] * jscale );
xtopleft = xcentre - sr;
ytopleft = ycentre - sr;
hr = 2 * sr;
// g.setColor( Color.darkGray );
// g.drawOval( ytopleft, xtopleft, hr, hr );
int xtoa = xtopleft;
int ytoa = ytopleft;
Color colour = Color.white;
Color lastColour = Color.white;
for ( int n=0; n< ap.column[ ap.nviewLatitude ].graphicLayerCount; n++ ) {
colour = colorpalette[ (int)ap.column[ ap.nviewLatitude ].graphicLayer[n][8] ];
// colour = setColour( (int)column[ nviewLatitude ].graphicLayer[n][MATERIALTYPE], column[ nviewLatitude ].graphicLayer[n][PHASE] );
if ( colour != lastColour ) {
r = ap.column[ ap.nviewLatitude ].graphicLayer[n][RADIUS];
if ( r > ap.column[ ap.nviewLatitude ].topOfAtmosphere ) r = ap.column[ ap.nviewLatitude ].topOfAtmosphere;
sr = (int)( r * jscale );
xtopleft = xcentre - sr;
ytopleft = ycentre - sr;
hr = 2 * sr;
g.setColor( colour );
g.fillOval( xtopleft, ytopleft, hr, hr );
lastColour = colour;
}
}
if ( showLayerBoundaries ) {
for ( int n=0; n<ap.column[ ap.nviewLatitude ].graphicLayerCount; n++ ) {
sr = (int)( ap.column[ ap.nviewLatitude ].graphicLayer[n][RADIUS] * jscale );
xtopleft = xcentre - sr;
ytopleft = ycentre - sr;
hr = 2 * sr;
g.setColor( Color.white );
g.drawOval( xtopleft, ytopleft, hr, hr );
}
}
// g.setColor( Color.darkGray );
// g.fillRect( 100, 120, 30, 30 );
paintTemperatureGraph( g, asteroid, ytoa, jscale );
}
void paintTemperatureGraph( Graphics g, GeologicalColumn asteroid, int xtoa, double jscale ) {
int count = 0;
int x1, x2, y1, y2, m;
double t, r, d;
ConductiveLayer layer;
ListIterator<ConductiveLayer> itr = null;
itr = ap.column[ ap.nviewLatitude ].llayer.listIterator();
x1 = 0;
y1 = 0;
t = 0;
g.setColor( Color.black );
// g.drawLine( centre-bradius, centre, centre, centre );
g.setColor( Color.yellow );
for ( int n=0; n<ap.column[ ap.nviewLatitude ].nlayers; n++ ) {
t = ap.column[ ap.nviewLatitude ].graphicLayer[n][4];
y1 = xtoa + (int)(( ap.column[ ap.nviewLatitude ].graphicLayer[0][2] - ap.column[ ap.nviewLatitude ].graphicLayer[n][2] ) * jscale );
x1 = (int)( ap.column[ ap.nviewLatitude ].graphicLayer[n][4] / 50.0 );
m = n+1;
if ( m == ap.column[ ap.nviewLatitude ].nlayers ) m = n;
y2 = xtoa + (int)(( ap.column[ ap.nviewLatitude ].graphicLayer[0][2] - ap.column[ ap.nviewLatitude ].graphicLayer[m][2] ) * jscale );
x2 = (int)( ap.column[ ap.nviewLatitude ].graphicLayer[m][4] / 50.0 );
g.drawLine( x1, y1, x2, y2 );
}
// g.drawLine( x1, y1, centre, y1 );
g.drawString( " " + t + " degrees K", x1, y1 );
g.drawString( " 0 degrees K", x1, (ymax/2) );
}
// paint 18 columbs
void paint( Graphics g, double topRockRadius ) {
double r;
int x, phase, ci, xview, xloc;
int xstep = (int)( (double)xmax / 18.0 ); // column pixel width
int xspace = xmax - 18 * xstep - 7; // left margin pixel width
boolean drawColSeparators = true;
boolean widenViewLatitude = false;
// show top layers
double topLayerDepth = 50000.0;
double jscale = xmax / topLayerDepth; // screen to show top 40 km
jscale = 20.0 * jscale;
// clear window
g.setColor( Color.white );
g.fillRect( 0, 0, xmax, ymax );
x = xmax-1;
// circle centre (pixels)
xcentre = xmax / 2;
ycentre = (int)( topRockRadius * jscale ) + ( ymax / 2 );
// paint atmosphere filled circle
// surface rock radius
int sr = (int)( ap.column[ 0 ].graphicLayer[0][RADIUS] * jscale );
for ( int nlat=0; nlat<ap.nLatitudes; nlat++ ) {
Color colour = Color.white;
Color lastColour = Color.white;
// draw vertical coloumn separator
if ( (nlat % 3) == 0 && drawColSeparators ) {
g.setColor( Color.lightGray );
g.drawLine( x, 0, x, ymax );
// if ( ap.year < 200 ) System.out.println( "nlat " + nlat + " x " + x + " xspace " + xspace );
x--;
colour = Color.lightGray;
}
x -= xstep;
for ( int n=ap.column[ nlat ].graphicLayerCount-1; n>=0; n-- ) {
if ( ap.column[ nlat ].layerTable[n].genericMaterialType() == 2 ) {
colour = setRockHeatColour( ap.column[ nlat ].layerTable[n].currentTemperatureK, ap.column[ nlat ].layerTable[n].initialTemperature, (int)20 );
} else {
phase = ap.column[ nlat ].layerTable[n].phase;
ci = (int)ap.column[ nlat ].layerTable[n].material.colourIndex[ phase ];
colour = colorpalette[ ci ];
}
// layer.material.colourIndex[ layer.phase ];
if ( colour != lastColour ) {
r = ap.column[ nlat ].layerTable[n].radius ;
if ( r > ap.column[ nlat ].topOfAtmosphere ) r = ap.column[ nlat ].topOfAtmosphere;
sr = (int)( r * jscale );
xtopleft = xcentre - sr;
ytopleft = ycentre - sr;
hr = 2 * sr;
// if ( ap.year < 100 && nlat == 9 )System.out.println( ap.year + " nBand " + ap.column[ nlat ].layerTable[n].nBand + " mtl " + ap.column[ nlat ].layerTable[n].material.index + " r " + (float)r + " ytop " + ytopleft + " sr " + sr );
// if ( ap.year < 100 && nlat == 9 )System.out.println( ap.year + " nBand " + ap.column[ nlat ].layerTable[n].nBand + " mtl " + ap.column[ nlat ].layerTable[n].material.index + " r " + (float)r + " radius " + ap.column[ nlat ].layerTable[n].radius + " ci " + ci );
// colour = colorpalette[ (int)( Math.random() * 5 ) ];
g.setColor( colour );
// if ( ap.column[ nlat ].layerTable[n].currentTemperatureK > 280 ) g.setColor( Color.red );
// g.fillOval( xtopleft, ytopleft, hr, hr );
xview = 0;
if ( ( nlat == ap.nviewLatitude ) && widenViewLatitude ) {
xview = xspace; // widen viewLatititude
x -= xview;
}
g.fillRect( x, ytopleft, xstep+xview, hr );
lastColour = colour;
}
}
}
}
// rock heat colour rises from brown to red and then yellow
// or it falls from brown to purple and then bluw
Color setRockHeatColour( double rockTemperatureK, double rockInitialTemperatureK, double scale ) {
Color heatColour = new Color( 90, 0, 0 );
double v = ( rockTemperatureK - rockInitialTemperatureK ) * scale;
int i;
int plus1 = 255 - 90; // red limit
int plus2 = plus1 + 255; // yellow limit
int minus1 = -90; // purple limit
int minus2 = minus1 - 90; // blue limit
if ( v > 0 ) {
if ( v < plus1 ) {
heatColour = new Color( 90 + (int)v, 0, 0 );
} else if ( v < plus2 ) {
i = (int)v - plus1;
if ( i >=0 && i <=255 ) {
heatColour = new Color( 255, i, 0 );
} else {
heatColour = new Color( 255, 127, 0 );
}
} else {
heatColour = new Color( 255, 255, 0 );
}
} else {
if ( v > minus1 ) {
heatColour = new Color( 90, 0, -(int)v );
} else if ( v > minus2 ) {
i = 90 + (int)v;
if ( i >=0 && i <= 255 ) {
heatColour = new Color( i, 0, 90 );
} else {
heatColour = new Color( 45, 0, 90 );
}
} else {
heatColour = new Color( 0, 0, 90 );
}
}
return( heatColour );
}
// rock heat colour rises from grey to red
// or it falls from grey to pale bluw
Color setRockHeatColour( double rockTemperatureK, double rockInitialTemperatureK, int range ) {
int gray = 90;
int r, g, b;
Color heatColour = new Color( 90, 0, 0 );
int t = (int)( rockTemperatureK - rockInitialTemperatureK );
if ( t > 0 ) {
if ( t > range ) t = range;
r = gray + ( ( 255 - gray ) * t ) / range;
g = gray - ( gray * t ) / range;
b = g;
} else {
if ( t < -range ) t = -range;
r = gray + ( gray * t ) / range;
g = gray - ( ( 255 - gray ) * t ) / range;
b = g;
}
if ( ap.year < 100 && t != 0 ) System.out.println( "T " + t + " RGB " + r + ", " + g + ", " + b );
heatColour = new Color( r, g, b );
return( heatColour );
}
// set colour using material index
Color setColour( int index, double phase ) {
Color colour = Color.white;
Material mt = new Material( ap.column[ ap.nviewLatitude ], index );
int cindex = 0;
if ( phase == 1 ) cindex = 1;
if ( phase == 2 ) cindex = 2;
colour = colorpalette[ mt.colourIndex[ cindex ] ];
return( colour );
}
}
// handle annual scheduled events (adapted from Orbit3D.java )
class EventHandler {
Glaciation ap;
long eventQueue[] = new long[500];
int nextEvent = 0;
int eventCount = 0;
int repeatEventType = 0;
long repeatEventCZYN;
long repeatEventTotal[] = new long[10];
int repeatEventCount;
long repeatEventStartCZYN;
long repeatEventStopCZYN;
int nrepeatEvent = 0;
long repeatEventMeasure[][] = new long[4][1000];
int repeatEventBody = 11;
long repeatEventInterval = 0; // quarter hour intervals of day
double timestepRemainder;
double secsperday;
long printAsteroidCZYN;
long setBandNumCZYN;
long rollingGraphCZYN;
long printEventQueueCZYN;
long layerDisplayCZYN;
long precipitstionCZYN;
long pauseEventCZYN;
long milankovitchCZYN;
// default constructor
public EventHandler() {}
public EventHandler( Glaciation a ) {
ap = a;
printAsteroidCZYN = (long)ap.startCZYN;
this.insertEvent( printAsteroidCZYN );
System.out.println( "printAsteroidCZYN " + (long)printAsteroidCZYN);
setBandNumCZYN = (long)ap.startCZYN + 9;
this.insertEvent( setBandNumCZYN );
System.out.println( "setBandNumCZYN " + (long)setBandNumCZYN);
layerDisplayCZYN = (long)ap.startCZYN + 10;
this.insertEvent( layerDisplayCZYN );
rollingGraphCZYN = (long)ap.startCZYN + 10;
this.insertEvent( rollingGraphCZYN );
precipitstionCZYN = (long)ap.startCZYN + ap.delayedStart; // delayed start of
// precipitstionCZYN = ap.eemianCZYN; // immediate start of precipitation
this.insertEvent( precipitstionCZYN );
System.out.println( ap.year + ": precipitationnCZYN " + precipitstionCZYN );
if ( ap.pauseDate > 0 ) {
this.pauseEventCZYN = (long)ap.startCZYN + (long)ap.pauseDate;
this.insertEvent( pauseEventCZYN );
}
if ( ap.solar == 3 ) {
this.milankovitchCZYN = (long)ap.startCZYN;
this.insertEvent( this.milankovitchCZYN );
}
System.out.println( "printAsteroidCZYN " + printAsteroidCZYN );
System.out.println( "setBandNumCZYN " + setBandNumCZYN );
System.out.println( "layerDisplayCZYN " + layerDisplayCZYN );
System.out.println( "rollingGraphCZYN " + rollingGraphCZYN );
System.out.println( "pauseEventCZYN " + pauseEventCZYN );
// System.out.println( "precipitstionCZYN " + precipitstionCZYN );
}
// a variety of events can all happen at the same time,
// so when an event occurs, a whole bunch of things need to be checked.
void checkEventQueue() {
double dt = 1.0;
double ta; // air temperature
double hl; // added layer height
double ri; // repeat interval
double hmax; // max snow/ice sheet depth
double lastDate;
boolean repaint = false;
if ( nextEvent < 500 ) {
// if ( ap.currentCZYN > 131885997 && ap.currentCZYN < 131886003 ) {
// System.out.println( ap.currentCZYN + " eventQueue[nextEvent] " + eventQueue[nextEvent] + " " + nextEvent + " :" + ap.year );
// printEventQueue();
// }
if ( (long)ap.currentCZYN == eventQueue[nextEvent] ) {
// if ( ap.printManager ) System.out.println( ap.currentCZYN + " EVENT at " + ap.currentCZYN );
if ( ap.printManager ) System.out.println( ap.currentCZYN + " EVENT at " + ap.currentCZYN );
if ( (long)ap.currentCZYN == pauseEventCZYN ) {
ap.go = false;
System.out.println( "paused at year "+ pauseEventCZYN );
// update display
ap.column[ ap.nviewLatitude ].setBandNumbers(); // get data for graphic display
ap.rollingGraphOutput( ap.century );
ap.rollingGraphUpdated = true;
repaint = true;
pauseEventCZYN += ap.pauseDate;
this.insertEvent( pauseEventCZYN );
}
if ( ap.currentCZYN == printAsteroidCZYN ) {
ap.column[ ap.nviewLatitude ].printhcmColumn( 0 );
if ( ap.year >= 11999 && ap.year <= 11998 ) {
printAsteroidCZYN += 1;
} else {
printAsteroidCZYN += 6000;
}
this.insertEvent( printAsteroidCZYN );
// if ( ap.printManager ) System.out.println( "print asteroid " );
// this.printEventQueue();
System.out.println( "next print asteroid date " + printAsteroidCZYN );
}
if ( ap.currentCZYN == setBandNumCZYN ) {
ap.column[ ap.nviewLatitude ].setBandNumbers();
setBandNumCZYN += 1000;
this.insertEvent( setBandNumCZYN );
if ( ap.printManager ) System.out.println( "set band numbers " );
}
/*
// if ( ( ap.currentCZYN == (ap.eemianCZYN) ) || ( ap.currentCZYN == (ap.eemianCZYN + 15018) ) ) {
if ( ap.currentCZYN == precipitstionCZYN ) {
// ap.column[ nviewLatitude ].addLayer( ap.column[ nviewLatitude ].thinSurfaceBand, ap.column[ nviewLatitude ].surfaceBandTemperature, 1000, 3 );
// end precipitation after 10000 years
// snow temperature is air temperature.
// if air temperature > 273K, snow temperature =
ta = ap.column[ ap.nviewLatitude ].airLowestBandTemperatureK;
if ( ta > 273 ) ta = 273;
hl = 25.0; // layer height
ri = 300.0; // deposition interval (years)
hmax = 500.0; // maximum height of layer
// if ( ap.currentCZYN < ap.eemianCZYN + 990010 ) {
// System.out.println( "CZYN " + (int)ap.currentCZYN + " " + (int)ap.eemianCZYN );
// ap.asteroid.addLayer( ap.asteroid.thinSurfaceBand, ap.asteroid.atmosphereTemperature, 200.0, 3 );
// ap.asteroid.addLayer( ap.asteroid.thinSurfaceBand, 260, 50, 3.1 );
// ap.asteroid.addLayer( ap.asteroid.thinSurfaceBand, 270, 500, 3.0 );
// ap.asteroid.addLayer( ap.asteroid.thinSurfaceBand, 270.00, 3000.0, 3.00 );
// ap.column[ ap.nviewLatitude ].addLayer( ap.column[ ap.nviewLatitude ].thinSurfaceBand, 270.00, 260.00, 3.05 );
// if ( ( ap.year < 50000 || ap.year > 100000 ) && ap.asteroid.iceWaterDepth < 250 ) {
// if ( ap.column[ ap.nviewLatitude ].iceWaterDepth < 250 && ap.year < 100000 ) {
if ( ap.column[ ap.nviewLatitude ].iceWaterDepth < hmax ) {
// System.out.println( "start adding glaciation in year " + ap.year );
// ap.column[ ap.nviewLatitude ].addIceLayer( ta, hl, 3.05 );
ap.column[ ap.nviewLatitude ].addLayer( ap.column[ ap.nviewLatitude ].thinSurfaceBand, ta, hl, 3.05 );
// ap.column[ ap.nviewLatitude ].addLayer( ap.column[ ap.nviewLatitude ].thinSurfaceBand, ta, hl, 12.0 );
// ap.column[ ap.nviewLatitude ].addIceLayer( 270.0, 10.00, 3.05 );
// ap.column[ ap.nviewLatitude ].addIceLayer( 270.0, 10.00, 3.05 );
// ap.column[ ap.nviewLatitude ].addIceLayer( 270.0, 10.00, 3.05 );
// ap.column[ ap.nviewLatitude ].addIceLayer( 270.0, 10.00, 3.05 );
// ap.column[ ap.nviewLatitude ].addIceLayer( 270.0, 10.00, 3.05 );
// System.out.println( ap.year + " layer height " + hl + " layer T " + ta );
// System.out.println();
}
ap.column[ ap.nviewLatitude ].startOfGlaciation = ap.year;
// ap.asteroid.setAirEmissivity = 0.999;
// ap.asteroid.flagSetAirEmissivity = true;
// ap.icedropCounter = 100;
// ap.solar = 0;
// precipitstionCZYN += 5000;
if ( ap.column[ ap.nviewLatitude ].varyIceAlbedo ) {
ap.column[ ap.nviewLatitude ].varyIceAlbedoReset = true;
}
ap.column[ ap.nviewLatitude ].resynchronise = true; // need to renumber layers
// } else {
// ap.column[ nviewLatitude ].addLayer( ap.column[ nviewLatitude ].thinSurfaceBand, ap.column[ nviewLatitude ].atmosphereTemperature, 25.0, 3 );
precipitstionCZYN += ri;
// }
this.insertEvent( precipitstionCZYN );
}
*/
// add snow to active columns
for ( int ncol=0; ncol<=17; ncol++ ) {
if ( ap.column[ ncol ].columnActiveCalculation ) {
if ( ap.column[ ncol ].snowStartCZYN <= ap.currentCZYN && ap.column[ ncol ].snowStopCZYN > ap.currentCZYN ) {
if ( ap.currentCZYN == ap.column[ ncol ].precipitstionCZYN ) {
System.out.println( ap.year + ": CZYN " + ap.currentCZYN );
ta = ap.column[ ncol ].airLowestBandTemperatureK; // air temperature
// if ( ta > 273 ) ta = 273;
hl = ap.column[ ncol ].snowLayerThickness; // layer height
ri = ap.column[ ncol ].snowDepositionInterval; // deposition interval (years)
hmax = ap.column[ ncol ].snowMaximumDepth; // maximum height of layer
if ( ap.column[ ncol ].iceWaterDepth < hmax ) {
ap.column[ ncol ].addLayer( ap.column[ ncol ].thinSurfaceBand, ta, hl, ap.column[ ncol ].snowType );
ap.column[ ncol ].resynchronise = true; // need to renumber layers
ap.column[ ncol ].precipitstionCZYN += ri;
this.insertEvent( (long)ap.column[ ncol ].precipitstionCZYN );
// if ( !ap.snowing ) ap.setInitialTemperatures();
ap.snowing = true;
} else {
ap.column[ ncol ].precipitstionCZYN += ri;
this.insertEvent( (long)ap.column[ ncol ].precipitstionCZYN );
}
}
}
}
}
if ( ap.currentCZYN == layerDisplayCZYN ) {
// if ( ap.currentCZYN - (long)ap.eemianCZYN > 65000000 ) System.out.println( ap.currentCZYN + " layerDisplay layerDisplayCZYN " + layerDisplayCZYN );
repaint = true;
layerDisplayCZYN += 100000;
this.insertEvent( layerDisplayCZYN );
}
if ( ap.currentCZYN == rollingGraphCZYN ) {
// if ( ap.albedoTest ) ap.column[ ap.nviewLatitude ].incrementThinSurefaceLayerAlbedo( );
ap.column[ ap.nviewLatitude ].setBandNumbers(); // get data for graphic display
ap.rollingGraphOutput( ap.century );
ap.rollingGraphUpdated = true;
repaint = true;
rollingGraphCZYN += ap.canvas1.calendarStep;
this.insertEvent( rollingGraphCZYN );
if ( ap.printManager ) System.out.println( (float)(ap.currentCZYN - ap.eemianCZYN) + " rollingGraph numbers " );
// if ( ap.currentCZYN > 131885700 && ap.currentCZYN < 131886003 ) System.out.println( ap.currentCZYN + " new rollingGraphCZYN " + rollingGraphCZYN );
}
if ( ap.currentCZYN == milankovitchCZYN ) {
//Milankovitch ebentManager returns next date
lastDate = ap.currentCZYN;
milankovitchCZYN = ap.mil.eventManager( milankovitchCZYN );
System.out.println( "Milankovitvh current " + ap.currentCZYN + " next " + milankovitchCZYN );
if ( lastDate !=milankovitchCZYN ) this.insertEvent( milankovitchCZYN ); // don't add event if same date
}
removeEvent( (long)ap.currentCZYN );
if ( repaint == true ) ap.repaint();
}
} else {
nextEvent = 0; // prevent overflow
}
}
void setupRepeatEvent( int type, double startCZYN, int bn, double duration, double interval ) {
repeatEventType = type;
repeatEventBody = bn;
if ( type == 1 ) {
repeatEventCZYN = (long)ap.currentCZYN;
repeatEventTotal[0] = 0;
repeatEventTotal[1] = 0;
repeatEventTotal[2] = 0;
repeatEventCount = 0;
repeatEventStopCZYN = repeatEventCZYN + (long)duration; // record for some period of time
repeatEventInterval = (long)interval;
repeatEventCZYN = ap.currentCZYN + (long)( interval ); // add quarter hour
insertEvent( repeatEventCZYN );
System.out.println( "currentCZYN " + ap.currentCZYN + " repeatEventCZYN scheduled for " + repeatEventCZYN );
printEventQueue();
} else if ( type == 2 || type == 3 || type == 4 || type == 5 ) {
repeatEventStartCZYN = (long)startCZYN;
repeatEventStopCZYN = (long)startCZYN + (long)duration;
repeatEventCount = 0;
repeatEventInterval = (long)interval;
repeatEventCZYN = (long)startCZYN;
insertEvent( (long)startCZYN );
System.out.println( "repeat event scheduled for " + startCZYN );
}
}
// A variety of dates need to be accurately hit. These include:
// 1. Historical dates on which known planetary body positions are corrected.
// 2. Stop dates when state vectors are printed.
// 3. Body state changes defined by b[n].initialCZYN and finalCZYN
// 4. Rockcloud halo launch dates
// The timestep manager keeps track of the nearest upcoming date from a variety of sources,
// and, if necessary, alters dt in order to hit these dates.
double timestepManager( double currentTimestep ) {
double nextjd, timestepToEvent;
double newTimestep = currentTimestep;
// if ( pool.currentCZYN > 2454877.2 && pool.currentCZYN < 2454877.4) System.out.println( "***MANAGER " + pool.currentCZYN + " nextCZYN : " + eventQueue[nextEvent] + " dt " + currentTimestep );
// System.out.println( "***MANAGER " + pool.currentCZYN + " nextCZYN : " + eventQueue[nextEvent] + " dt " + currentTimestep );
// if there is a remaining timestep to complete, use it as the next timestep
// so as to not break the timestep sequence.
if ( timestepRemainder != 0.0 ) newTimestep = timestepRemainder;
nextjd = eventQueue[nextEvent];
// if ( nextjd <= pool.currentCZYN ) nextEvent++; // if next event is before current date, try next event along.
// if timestep length < current dt, return shortest.
timestepToEvent = ( nextjd - ap.currentCZYN ) * secsperday;
if ( Math.abs(timestepToEvent) < Math.abs(newTimestep) ) {
timestepRemainder = newTimestep - timestepToEvent;
newTimestep = timestepToEvent;
// if ( printManager ) {
// if ( newTimestep != dt ) System.out.println( "MANAGER " + pool.currentCZYN + " new timestep : " + newTimestep + " current timestep : " + currentTimestep + " rem " + timestepRemainder );
// }
} else {
timestepRemainder = 0;
}
// if ( pool.currentCZYN > 2454877.2 && pool.currentCZYN < 2454877.4) System.out.println( "MANAGER " + pool.currentCZYN + " nextjd " + nextjd + " timestep : " + newTimestep + " rem " + timestepRemainder );
return( newTimestep );
}
double getNearestDate( double currentDate, double currentNearestDate, double offeredDate, double timestep ) {
double newDate = currentNearestDate;
if ( timestep > 0 ) {
if ( offeredDate > currentDate ) {
if ( offeredDate < currentNearestDate ) newDate = offeredDate;
}
} else {
if ( offeredDate < currentDate ) {
if ( offeredDate > currentNearestDate ) newDate = offeredDate;
}
}
return( newDate );
}
// add event to event queue if it doesn't already exist
void addEvent( long jd ) {
int n;
boolean found = false;
for ( n=0; n<eventCount; n++ ) {
if ( eventQueue[n] == jd ) {
found = true;
break;
}
}
if ( !found ) {
eventQueue[ eventCount ] = jd;
eventCount++;
}
}
// Insert Julian date jd into date-ordered event queue
void insertEvent( long jd ) {
int n, i;
boolean found = false;
if ( eventCount == 0 ) {
eventQueue[0] = jd; // if no events, add new event
eventCount++;
// eventQueue[1] = jd * 2; // dummy terminal event
// eventCount++;
} else {
for ( n=0; n<eventCount; n++ ) {
if ( eventQueue[n] == jd ) {
found = true; // event already exists
break;
} else if ( eventQueue[n] > jd ) {
eventCount++;
for ( i=eventCount; i>n; i-- ) {
eventQueue[i] = eventQueue[i-1]; // move subsequent events down queue
}
eventQueue[n] = jd; // add new event
break;
} else if ( n == eventCount-1 ) {
eventCount++; // add new event to end of the queue
eventQueue[n+1] = jd;
}
}
}
// if ( ap.currentCZYN > 131885000 && ap.currentCZYN < 131886003 ) {
// System.out.println( ap.currentCZYN + " new event inserted at " + jd );
// printEventQueue();
// }
}
// Insert Julian date jd into date-ordered event queue before NextEvent pointer
void insertNowEvent( long jd ) {
for ( int i=eventCount; i>nextEvent; i-- ) {
eventQueue[i] = eventQueue[i-1]; // move subsequent events down queue
}
eventQueue[nextEvent] = (long)jd; // add new event
eventCount++;
nextEvent++; // point to next event after it
}
// Insert Julian date jd into date-ordered event queue
void insertIntermediateEvent( long jd ) {
int n, i;
boolean found = false;
if ( jd > eventQueue[0] && jd < eventQueue[eventCount-1]) {
for ( n=0; n<eventCount; n++ ) {
if ( eventQueue[n] == jd ) {
found = true; // event already exists
break;
} else if ( eventQueue[n] > jd ) {
eventCount++;
for ( i=eventCount; i>n; i-- ) {
eventQueue[i] = eventQueue[i-1]; // move subsequent events down queue
}
eventQueue[n] = (long)jd; // add new event
break;
}
}
}
}
// remove Julian date jd event from event queue
void removeEvent( long jd ) {
int n, i;
boolean found = false;
for ( n=0; n<eventCount; n++ ) {
if ( eventQueue[n] == jd ) {
for ( i=n; i<eventCount; i++ ) {
eventQueue[i] = eventQueue[i+1]; // move subsequent events up queue
eventQueue[i+1] = 0;
}
eventCount--;
break;
}
}
}
void printEventQueue() {
System.out.println( "Event Queue count = " + eventCount );
for ( int n=0; n<eventCount; n++ ) {
System.out.println( n + " " + eventQueue[n] );
}
System.out.println( "nextEvent " + nextEvent );
}
}
class SolarRadiation {
// source: Solar Radiation Calculation Dr. Mohamad Kharseh
// source: https://www.researchgate.net solarRadiationCalc.pdf by Dr. Mohamad Kharseh
// mean terrestrial solar radiation of 342 W/m^2 = 2.95E7 W / m^2 day
// Earth Mean Orbital Elements (J2000) source: https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html
// Semimajor axis (AU) 1.00000011
// Orbital eccentricity 0.01671022
// Orbital inclination (deg) 0.00005
// Longitude of ascending node (deg) -11.26064
// Longitude of perihelion (deg) 102.94719
// Mean Longitude (deg) 100.46435
double timeStepTable[][] = new double[2][370]; // J
int nTimeSteps;
double timeStep; // days
double meanSolarPower;
double latitude; // +/-degrees
double solarConstant; // Watts / square metre
int currentTimeStep;
double eccentricity;
double obliquity; // radians
double perihelionLongitude; // radians
double secondsPerYear = 3.154e+7;
public SolarRadiation() { }
// test code
public SolarRadiation( double solarConstant ) {
timeStepTerrestrialRadiation( solarConstant, 183.75, 185.25, 65 );
}
// create table of solar heat gains during each time step at given latiude
// try to set integer number of timesteps per year (e.g. 1000 per year)
public SolarRadiation( Glaciation a, GeologicalColumn gc, double solarConstant, double latitude, double timestepYears, double daysPerYear, double e, double ob, double pl ) {
Glaciation ap = a;
double ho;
double timestepDays = daysPerYear * timestepYears;
int nTimestepsPerYear = (int)( daysPerYear / timestepDays );
this.solarConstant = solarConstant;
this.latitude = latitude;
this.timeStep = timestepDays;
this.nTimeSteps = nTimestepsPerYear;
this.eccentricity = e;
this.obliquity = ob;
this.perihelionLongitude = pl;
this.currentTimeStep = 0;
double daynum = 0;
double totho = 0;
int step = 1;
timeStepTable[0][ 0 ] = 0;
for ( int s=1; s<=nTimestepsPerYear; s++ ) {
ho = timeStepTerrestrialRadiation( solarConstant, daynum, daynum+timestepDays, latitude );
totho += ho;
timeStepTable[1][ step ] = totho;
daynum += timestepDays;
timeStepTable[0][ step ] = ho;
step++;
}
this.meanSolarPower = totho / ap.secondsPerYear; // watts
// System.out.println( "total solar gain " + (long)totho + " mean power (watts) " + this.meanSolarPower );
System.out.println( "nTimestepsPerYear " + nTimestepsPerYear );
timeStepTable[0][ step ] = -1; // terminator character if needed
// printTimeStepTable();
gc.meanTimestepSolarGain = totho / (double)nTimestepsPerYear;
// System.out.println(" totho " + totho + " steps " + nTimestepsPerYear + " meanTimestepSolarGain " + (int)gc.meanTimestepSolarGain );
this.meanSolarPower = totho / this.secondsPerYear; // assuming time step one day
// System.out.println(" mean annual solar power " + this.meanSolarPower );
}
void printTimeStepTable() {
int n = 0;
System.out.println( "timeStepTable at latitude " + this.latitude );
while ( this.timeStepTable[0][n] >= 0 ) {
System.out.printf("% 9d", n );
System.out.printf("% 9d", (int)this.timeStepTable[0][n] );
System.out.println();
n++;
}
}
// return the solar gain (J) for the current time step, and update the current time step number
double getNextSolarGain() {
double solargain = this.timeStepTable[0][ this.currentTimeStep ];
this.currentTimeStep++;
if ( this.timeStepTable[0][ this.currentTimeStep ] < 0 ) this.currentTimeStep = 0;
return( solargain );
}
// return total solar gain received between startDay and endDay
// endDay > startDay in range 0 to 365
double getTotalSolarGain( int startDay, int endDay ) {
// System.out.println( "startDay" + startDay + " endDay " + endDay );
double solargain = this.timeStepTable[1][endDay] - this.timeStepTable[1][startDay];
return( solargain );
}
// tested, working for non-integer timesteps of n days.
// use to create a timestepTable of nsteps of terrestrial solar heat gains at different latitudes
double timeStepTerrestrialRadiation( double solarConstant, double dn1, double dn2, double latitude ) {
double toth = 0;
double h = 0;
// parse start and stop day numbers
if ( dn2 > dn1 ) {
// start day
boolean s1i = false;
int sd1 = (int)dn1;
double f11 = dn1 - (int)dn1; // remaining fraction of first day
if ( f11 == 0 ) s1i = true; // integer start day
f11 = fractionAngle( f11 );
double f12 = 180;
// stop day
boolean s2i = false;
int sd2 = (int)dn2;
double f22 = dn2 - (int)dn2;
if ( f22 == 0 ) s2i = true; // integer stop day
f22 = fractionAngle( f22 ); // initial fraction of last day
double f21 = -180;
// System.out.println( dn1 + " to " + dn2 + ": " + sd1 + " " + s1i + " + " + sd2 + " " + s2i );
if ( sd1 == sd2 ) {
// period is fraction of one day, sd1
f12 = f22;
// System.out.println( sd1 + ": " + f11 + " to " + f12 + " same day" );
h = hourlyExtraTerrestrialRadiation( solarConstant, sd1, f11, f12, latitude );
// System.out.println( "h = " + (int)h );
toth += h;
} else {
int wd1 = sd1 + 1;
int wd2 = sd2 - 1;
if ( !s1i ) {
// start day fraction
// System.out.println( sd1 + ": " + f11 + " to " + f12 + " start" );
h = hourlyExtraTerrestrialRadiation( solarConstant, sd1, f11, f12, latitude );
// System.out.println( "h = " + (int)h );
toth += h;
} else {
wd1 = sd1;
}
// whole days
for ( int n=wd1; n<=wd2; n++ ) {
// System.out.println( n + " whole" );
h = dailyExtraTerrestrialRadiation( solarConstant, wd1, latitude );
// System.out.println( "h = " + (int)h );
toth += h;
}
// stop day fraction
if ( !s2i ) {
// System.out.println( sd2 + ": " + f21 + " to " + f22 + " stop" );
h = hourlyExtraTerrestrialRadiation( solarConstant, sd2, f21, f22, latitude );
// System.out.println( "h = " + (int)h );
toth += h;
} else {
wd1 = sd1;
}
}
// System.out.println( "toth = " + (int)toth );
} else {
System.out.println( "invalid day numbers");
}
return( toth );
}
// when fraction is 0.5, this corresponds to noon, hour angle 0
// when fraction is 0 or 1, hour angle = 180
double fractionAngle( double fraction ) {
double angle = ( fraction - 0.5 ) * 360.0;
// System.out.println( fraction + " a " + angle );
return( angle );
}
// degrees are coverted to radians
// tested working (produces slighly smaller values over a day than dailETradiation))
double hourlyExtraTerrestrialRadiation( double solarConstant, double dayNumber, double ha1, double ha2, double latitude ) {
double hourlyRadiation = 0;
double hourAngle1 = degreesToRadians( ha1 );
double hourAngle2 = degreesToRadians( ha2 );
double lat = degreesToRadians( latitude ); // radians
double declination = declination( dayNumber ); // radians
double sunsetHourAngle = hourAngle( lat, declination ); // radians
double hoursSunlight = 2.0 * radiansToDegrees( sunsetHourAngle ) / 15.0;
double f1 = 12.0 * 3600.0 * solarConstant / Math.PI;
// 0,033 is twice the eccentricity of the Earth's orbit
double f2 = 1.0 + (2.0 * this.eccentricity) * Math.cos ( degreesToRadians( ( 360.0 * dayNumber ) / 365.0 ) );
double f3 = Math.cos(lat) * Math.cos(declination)* ( Math.sin(hourAngle2) - Math.sin(hourAngle1 ) );
double f4 = ( hourAngle2 - hourAngle1 ) * Math.sin(lat) * Math.sin(declination);
hourlyRadiation = f1 * f2 * ( f3 + f4 );
if ( hourlyRadiation < 0 ) hourlyRadiation = 0;
return( hourlyRadiation );
}
// degrees are coverted to radians
// tested working
double dailyExtraTerrestrialRadiation( double solarConstant, double dayNumber, double latitude ) {
double dailyRadiation;
double lat = degreesToRadians( latitude ); // radians
double declination = declination( dayNumber ); // radians
double sunsetHourAngle = hourAngle( lat, declination ); // radians
double f1 = 24.0 * 3600.0 * solarConstant / Math.PI;
// 0,033 is twice the eccentricity of the Earth's orbit
double f2 = 1.0 + (2.0 * this.eccentricity) * Math.cos ( degreesToRadians( ( 360.0 * dayNumber ) / 365.0 ) );
double f3 = Math.cos(lat) * Math.cos(declination)* Math.sin(sunsetHourAngle );
double f4 = sunsetHourAngle * Math.sin(lat) * Math.sin(declination);
// was f4 = Math.PI * sunsetHourAngle(degrees) * Math.sin(lat) * Math.sin(declination ) / 180.0;
dailyRadiation = f1 * f2 * ( f3 + f4 );
if ( dailyRadiation <= 0 ) {
dailyRadiation = 0;
// System.out.println( latitude + " sunsetHourAngle " + sunsetHourAngle + " declination " + declination);
}
return( dailyRadiation );
}
double instantaneousExtraTerrestrialRadiation( double solarConstant, double dayNumber, double hourAngle, double latitude ) {
double instantaneousRadiation;
double hourAngleR = degreesToRadians( hourAngle );
double lat = degreesToRadians( latitude ); // radians
double declination = declination( dayNumber ); // radians
double f1 = solarConstant;
// 0,033 is twice the eccentricity of the Earth's orbit
double f2 = 1.0 + (2.0 * this.eccentricity) * Math.cos ( degreesToRadians( ( 360.0 * dayNumber ) / 365.0 ) );
double f3 = Math.cos(lat) * Math.cos(declination)* Math.cos(hourAngleR);
double f4 = Math.sin(lat) * Math.sin(declination);
instantaneousRadiation = f1 * f2 * ( f3 + f4 );
if ( instantaneousRadiation < 0 ) instantaneousRadiation = 0;
return( instantaneousRadiation );
}
// Declination is the angle made between the plane of the equator
// and the line joining the two centres of the earth and the sun
// tested working
double declination( double dn ) {
double ob = radiansToDegrees( this.obliquity );
// double d = ob * Math.sin( degreesToRadians( 360 * ( 284.0 + dn ) / 365.0 ) );
// dayOffset = 284 when this.perihelionLongitude = 1.7962.
// and has a range of 365 days when when this.perihelionLongitude is between 0 and 2.PI radians
double dayOffset = 284.0 + ( this.perihelionLongitude - 1.7962 ) * 365.25 / ( 2.0 *Math.PI );
double d = ob * Math.sin( degreesToRadians( 360 * ( dayOffset + dn ) / 365.0 ) );
return( degreesToRadians( d ) );
}
// The hour angle is the sun’s angular deviation from south
// tested working
double hourAngle( double lat, double dec ) {
double ha;
double rightAngle = 0.5 * Math.PI;
if ( lat >= 0 ) {
// northern hemisphere (tested)
if ( (lat-dec) > rightAngle ) {
// sun never rises
ha = 0;
} else if ( (lat+dec) > rightAngle ) {
// sun never sets
ha = Math.PI;
} else {
// sun rises and sets
ha = Math.acos( -Math.tan(lat) * Math.tan(dec) );
}
} else {
// southern hemisphere (untested)
if ( (-lat-dec) > rightAngle ) {
// sun never sets
ha = Math.PI;
} else if ( (-lat+dec) > rightAngle ) {
// sun never rises
ha = 0;
} else {
// sun rises and sets
ha = Math.acos( -Math.tan(lat) * Math.tan(dec) );
}
}
// if ( lat == rightAngle ) System.out.println( "ha " + ha );
return( ha );
}
double degreesToRadians( double d ) {
double r = d * Math.PI / 180.0;
return( r );
}
double radiansToDegrees( double r ) {
double d = r * 180.0 / Math.PI;
return( d );
}
// untested
double arcos(double x) {
double a = Math.atan( ( Math.sqrt(1-x*x) ) / x );
return( a );
}
private void printDebugData(JTable table, double step1 ) {
int numRows = table.getRowCount();
int numCols = table.getColumnCount();
javax.swing.table.TableModel model = table.getModel();
/*
source: https://www.cs.colostate.edu/~cs160/.Summer16/resources/Java_printf_method_quick_reference.pdf
Java printf( ) Method Quick Reference
System.out.printf( “format-string” [, arg1, arg2, … ] );
Format String:
Composed of literals and format specifiers. Arguments are required only if there are format specifiers in the
format string. Format specifiers include: flags, width, precision, and conversion characters in the following
sequence:
% [flags] [width] [.precision] conversion-character ( square brackets denote optional parameters )
Flags:
- : left-justify ( default is to right-justify )
+ : output a plus ( + ) or minus ( - ) sign for a numerical value
0 : forces numerical values to be zero-padded ( default is blank padding )
, : comma grouping separator (for numbers > 1000)
: space will display a minus sign if the number is negative or a space if it is positive
Width:
Specifies the field width for outputting the argument and represents the minimum number of characters to
be written to the output. Include space for expected commas and a decimal point in the determination of
the width for numerical values.
Precision:
Used to restrict the output depending on the conversion. It specifies the number of digits of precision when
outputting floating-point values or the length of a substring to extract from a String. Numbers are rounded
to the specified precision.
Conversion-Characters:
d : decimal integer [byte, short, int, long]
f : floating-point number [float, double]
c : character Capital C will uppercase the letter
s : String Capital S will uppercase all the letters in the string
h : hashcode A hashcode is like an address. This is useful for printing a reference
n : newline Platform specific newline character- use %n instead of \n for greater compatibility
Examples:
System.out.printf("Total is: $%,.2f%n", dblTotal);
System.out.printf("Total: %-10.2f: ", dblTotal);
System.out.printf("% 4d", intValue);
System.out.printf("%20.10s\n", stringVal);
String s = "Hello World";
System.out.printf("The String object %s is at hash code %h%n", s, s);
String class format( ) method:
You can build a formatted String and assign it to a variable using the static format method in the String class.
The use of a format string and argument list is identical to its use in the printf method. The format method
returns a reference to a String. Example:
String grandTotal = String.format("Grand Total: %,.2f", dblTotal);
*/
System.out.println("Value of data: ");
for (int i=0; i < numRows; i++) {
System.out.printf( "% 6.1f", 2.5 + (double)(i*step1) );
for (int j=0; j < numCols; j++) {
System.out.printf("% 9d", model.getValueAt(i, j) );
// System.out.print(" " + model.getValueAt(i, j));
}
System.out.println();
}
System.out.println("--------------------------");
}
}
class Milankovitch {
Glaciation exp;
long currentEpoch;
long irregularStartDate;
long irregularStopDate;
int irregularRowCount;
long regularInterval;
double eventTable[][] = new double[1000][4];
int eventTablerowcount;
double currentDate;
double eccentricity;
double obliquity;
double perihelionLongitude;
double nextDate;
double degreesToRadians = Math.PI / 180;
public Milankovitch() {
}
public Milankovitch( Glaciation a, String filename, long epoch, int interval ) {
this.exp = a;
this.currentEpoch = epoch;
this.regularInterval = interval;
System.out.println( "Milankovitch start" );
creatMergedTable( filename );
System.out.println( "Milankovitch end" );
}
void creatMergedTable( String filename ) {
double irregular[][] = csvReader( filename );
int index = 0;
int irregularIndex = 0;
int regularIndex = 0;
long irregularDate; // current regular date
long currentDate = irregularStartDate; // current date
long regularDate = irregularStartDate; // current irregular date
// first entry is first irrgular date
eventTable[index][0] = irregularStartDate;
eventTable[index][1] = irregular[irregularIndex][1];
eventTable[index][2] = irregular[irregularIndex][2];
eventTable[index][3] = irregular[irregularIndex][3];
currentDate = irregularStartDate;
System.out.println( "[" + index + "] " + (long)eventTable[index][0] + ": " + eventTable[index][1] + " " + eventTable[index][2] * this.degreesToRadians + " " + eventTable[index][3] * this.degreesToRadians );
// System.out.println( index + " merge idate " + irregularStartDate );
index++;
irregularIndex++;
regularDate += this.regularInterval;
while ( regularDate <= irregularStopDate ) {
irregularDate = (long)irregular[ irregularIndex ][0] + this.currentEpoch;
while ( regularDate >= irregularDate ) {
eventTable[index][0] = irregularDate;
eventTable[index][1] = irregular[irregularIndex][1];
eventTable[index][2] = irregular[irregularIndex][2];
eventTable[index][3] = irregular[irregularIndex][3];
currentDate = irregularDate;
System.out.println( "[" + index + "] " + (long)eventTable[index][0] + ": " + eventTable[index][1] + " " + eventTable[index][2] * this.degreesToRadians + " " + eventTable[index][3] * this.degreesToRadians );
// System.out.println( index + " merge idate " + irregularDate + " " + irregularIndex + " " + irregular[ irregularIndex ][0] );
index++;
irregularIndex++;
if ( irregularIndex > irregularRowCount ) break;
irregularDate = (long)irregular[ irregularIndex ][0] + this.currentEpoch;
}
if ( irregularIndex > irregularRowCount ) break;
// add emptry regular event to eventTable[]
if ( regularDate >= currentDate ) {
eventTable[index][0] = regularDate;
eventTable[index][1] = -1;
eventTable[index][2] = -1;
eventTable[index][3] = -1;
currentDate = regularDate;
System.out.println( "[" + index + "] " + (long)eventTable[index][0] + ": " + eventTable[index][1] + " " + eventTable[index][2] * this.degreesToRadians + " " + eventTable[index][3] * this.degreesToRadians );
// System.out.println( index + " merge date " + regularDate );
}
regularDate += this.regularInterval;
index++;
}
// System.out.println( "print 0-35" );
// for ( index=0; index<35; index++ ) {
// System.out.println( "[" + index + "] " + (long)eventTable[index][0] + ": " + eventTable[index][1] + " " + eventTable[index][2] + " " + eventTable[index][3] );
// }
// interpolate values to fill table
double v1, v2;
double t1, t2, t3;
int i[][] = {
{ 0, 0, 0, 0, },
{ 0, 0, 0, 0, }
};
index = 0;
while ( eventTable[ index ][0] < irregularStopDate ) {
for ( int j=1; j<=3; j++ ) {
if ( index == i[1][j] ) {
// re-index
i[0][j] = index;
i[1][j]++;
while ( eventTable[ i[1][j] ][j] == -1 ) {
i[1][j]++;
}
if ( eventTable[ i[0][j] ][0] < 64816000 ) System.out.println( j + ": " + i[0][j] + " " + i[1][j] + " v: " + eventTable[ i[0][j] ][j] + " " + eventTable[ i[1][j] ][j] );
}
v1 = eventTable[ i[0][j] ] [j];
v2 = eventTable[ i[1][j] ] [j];
t1 = eventTable[ i[0][j] ] [0];
t2 = eventTable[ i[1][j] ] [0];
t3 = eventTable[ index ] [0];
eventTable[ index ][j] = interpolate ( v1, t1, v2, t2, t3 );
}
if ( eventTable[ index ][0] < 67816000 ) System.out.println( "[" + index + "] " + (long)eventTable[ index ][0] + " " + (float)eventTable[ index ][1] + " " + (float)eventTable[ index ][2] * this.degreesToRadians + " " + (float)eventTable[ index ][3] * this.degreesToRadians );
index++;
}
eventTablerowcount = index;
/*
eventManager( 65938000 );
eventManager( 65938020 );
eventManager( 65693000 );
eventManager( 66111000 );
*/
}
double interpolate ( double v1, double t1, double v2, double t2, double t3 ) {
double v = v1;
double f;
if ( t2 > t1 ) {
f = ( t3 - t1 ) / ( t2 - t1 );
v = v1 + f * ( v2 - v1 );
}
return( v );
}
double[][] csvReader( String file ) {
String dataRow;
String rs[] = new String[10];
int nlayers, rowcount, fieldcount;
boolean bodyok;
boolean filefound = false;
BufferedReader CSVFile;
// double variable[] = new double[10];
double variable[][] = new double[500][4];
rowcount = 0;
try {
CSVFile = new BufferedReader(new FileReader(file));
rowcount = 0;
try {
dataRow = CSVFile.readLine(); // Read first line.
System.out.println( dataRow );
} catch (IOException e) {
dataRow = new String("");
}
while ( dataRow != null ) {
bodyok = false;
filefound = true;
StringTokenizer st = new StringTokenizer(dataRow, ",");
fieldcount = 0;
while (st.hasMoreTokens()) {
String s = st.nextToken();
rs[ fieldcount] = s;
fieldcount++;
}
if (rowcount > 0) {
variable[rowcount-1][0] = Double.valueOf( rs[0] ).doubleValue();
variable[rowcount-1][1] = Double.valueOf( rs[1] ).doubleValue();
variable[rowcount-1][2] = Double.valueOf( rs[2] ).doubleValue();
variable[rowcount-1][3] = Double.valueOf( rs[3] ).doubleValue();
if ( rowcount == 1 ) this.irregularStartDate = this.currentEpoch + (long)variable[rowcount-1][0];
System.out.println( rowcount + ": " + variable[rowcount-1][0] + ", " + variable[rowcount-1][1] + ", " + variable[rowcount-1][2] + ", " + variable[rowcount-1][3] + ", " );
}
rowcount++;
try {
dataRow = CSVFile.readLine(); // Read next line.
} catch (IOException e) {
dataRow = new String("");
}
}
// Close the file once all data has been read.
try {
CSVFile.close();
} catch (IOException e) {
}
} catch (IOException e) {
filefound = false;
System.out.println( file + " CSV file not found");
}
if ( filefound ) {
// End the printout with a blank line.
System.out.println( "rowcount " + rowcount + " variable[rowcount-1][0] " + (long)variable[rowcount-1][0] );
this.irregularStopDate = this.currentEpoch + (long)variable[rowcount-1][0];
this.irregularStopDate = 66098000; // correct value forced
System.out.println( "rowcount " + rowcount + " variable[rowcount-1][0] " + variable[rowcount-1][0] );
this.irregularRowCount = rowcount - 2;
System.out.println( this.currentEpoch + " " + variable[0][0] + " " + variable[rowcount-1][0] );
System.out.println("this.irregularStartDate " + this.irregularStartDate );
System.out.println("this.irregularStopDate " + this.irregularStopDate );
System.out.println("CSV File complete + " + rowcount + " rows " );
}
return( variable );
}
// get eccentricity, obliquity, perihelion longitude for cureent date
// if within date range sets current date's eccentricity, obliquity, and perihelion longitude
// else sets eccentricity, obliquity, and perihelion longitude to default values
long eventManager( long date ) {
double lat;
boolean dateFound = false;
GeologicalColumn gc;
double SRtimestep = 1.0 / 365.25;
double degreesToRadians = Math.PI / 180;
this.nextDate = this.currentDate;
// lookup date in eventTable
int index = 0;
if ( date >= irregularStartDate && date <= irregularStopDate ) {
while ( index < eventTablerowcount ) {
if ( eventTable[ index ][0] == date ) {
this.currentDate = eventTable[ index ][0];
this.eccentricity = eventTable[ index ][1];
this.obliquity = eventTable[ index ][2];
this.perihelionLongitude = eventTable[ index ][3] ;
if ( this.currentDate < this.irregularStopDate ) {
this.nextDate = eventTable[ index+1 ][0];
} else {
this.nextDate = this.irregularStopDate + 1;
}
dateFound = true;
} else if ( eventTable[ index ][0] > date ) {
// get previous index values
if ( index > 0 ) index = index - 1;
this.currentDate = eventTable[ index ][0];
this.eccentricity = eventTable[ index ][1];
this.obliquity = eventTable[ index ][2];
this.perihelionLongitude = eventTable[ index ][3] ;
if ( this.currentDate < this.irregularStopDate ) {
this.nextDate = eventTable[ index+1 ][0];
} else {
this.nextDate = this.irregularStopDate + 1;
}
dateFound = true;
}
if ( nextDate == date ) nextDate = eventTable[ index+2 ][0]; // fix for 2 dates the same
if ( dateFound ) break;
index++;
}
}
//
if ( date < this.irregularStartDate ) {
this.nextDate = eventTable[0][0]; //start of table
} else if ( date > this.irregularStopDate ) {
this.nextDate = this.irregularStopDate + 1; // off end of table
}
// get solar radiation data for all active latitudes
for ( int nlat=0; nlat<18; nlat++ ) {
lat = (nlat * 5) + 2.5;
gc = exp.column[ nlat ];
if ( dateFound ) {
// use Milankovitch data
if ( gc.columnActiveCalculation ) {
// use Milankovitch data
gc.sr = new SolarRadiation( this.exp, this.exp.column[ nlat ], 1367.0, lat, SRtimestep, 365.25, this.eccentricity, this.obliquity * degreesToRadians, this.perihelionLongitude * degreesToRadians );
}
} else {
if ( gc.columnActiveCalculation ) {
// use default data
gc.sr = new SolarRadiation( this.exp, this.exp.column[ nlat ], 1367.0, lat, SRtimestep, 365.25, 0.016702, 0.409093, 1.796257 );
}
}
}
if ( !dateFound ) {
System.out.println( "Milankovitch date " + date + " not found" );
this.nextDate = date; // return current date to flag end of eventTable
} else {
System.out.println( "Milankovitch date " + date + " nearest date is " + (long)this.currentDate );
}
return( (long)this.nextDate );
}
}
class Precipitation {
double meanGlobalPrecipitation = 0.814; // annua; precipiation, metres source http://www.mdpi.com/2073-4433/9/4/138/htm
double maxElevation = 3000; // max eelevation of rainfall, metres
double slope;
double latitude;
double asteroidRadius;
double degreesToRadians = Math.PI / 180.0;;
public Precipitation() {
}
public Precipitation( double lat, double maxelev ) {
latitude = lat;
maxElevation = maxelev;
slope = -meanGlobalPrecipitation / maxElevation;
}
double latitudeAdjustedPrecipitation( double elevation ) {
double pmean = annualMeanPrecipitation( elevation );
pmean = pmean * evaporationFactor( latitude );
return( pmean );
}
// linear decrease in precipiation with elevation
// y = mx + c where y = precipitation, x = elevation, c = meanGlobalPrecipitation
double annualMeanPrecipitation( double elevation ) {
double ap = slope * elevation + meanGlobalPrecipitation;
if ( ap < 0 ) ap = 0;
return( ap );
}
// Assuming all higher latitudes are glaciated, and all lower latitudes are unglaciated,
// water evaporation factor is unglaciated area/total araea
// area of arc of sphere = 2.pi.R.h
double evaporationFactor( double latitude) {
double lradians = latitude * degreesToRadians;
double e = Math.sin( lradians );
return( e );
}
}
class MaterialPhase {
// first find which side of vapour line the t-p point is.
// if vapour, return state 3
// else find which side of solid-liquid line t-p point is,
// if liquid, return state 2
// else state 3
// inputs layer temperature tK, layer pressure pPa
// returns state 0 = solid, 1-liquid, 2=vapour
// tested working 6 june 2018
public static double phase( double tK, double pPa ) {
double[][] water = {
// 0 1 2
// n, T, p,
{ 1, 0, 0 }, // minimum line
{ 1, 210, 1 }, // vapour line
{ 1, 230, 10 }, // vapour line
{ 1, 250, 100 }, // vapour line
{ 1, 273.15, 611.567 }, // triple point
{ 1, 300, 3500 }, // vapour line
{ 1, 320, 10e3 }, // vapour line
{ 1, 350, 35e3 }, // vapour line
{ 1, 373.35, 101.325e3 }, // boiling point at 1 atmosphere pressure
{ 1, 450, 0.95e6 }, // vapour line
{ 1, 550, 6.50e6 }, // vapour line
{ 1, 647, 22.064e6 }, // critical point
{ 1, 650, 10.00e12 }, // maximum line
{ 2, 273.16, 0 }, // minimum line
{ 2, 273.16, 611.567 }, // triple point
{ 2, 273.15, 101.325e3 }, // freezing point at 1 atmosphere pressure
{ 2, 251.165, 209.9e6 }, // solid-liquid line
{ 2, 256.164, 350.1e6 }, // solid-liquid line
{ 2, 272.990, 632.4e6 }, // solid-liquid line
{ 2, 355.000, 2.219e9 }, // solid-liquid line
{ 2, 650.000, 10.5e9 }, // solid-liquid line
{ 2, 650.000, 10.00e12 }, // maximum line
{ -1, }
};
double xmin = 1e20;
double xmax = -1e20;
double ymin = 1e20;
double ymax = -1e20;
double state = 0;
double nphaseline = 1; // vapour line is nphaseline 1
int n = 0;
Point3d p = new Point3d( tK, pPa, 0 );
Point3d p1 = new Point3d( water[n][1], water[n][2], 0);
Point3d p2 = new Point3d( 0, 0, 0);
Point3d v = new Point3d( 0, 0, 0 );
// findphase diagram max and min values
while ( water[n][0] > 0 ) {
if ( water[n][1] < xmin ) xmin = water[n][1];
if ( water[n][1] > xmax ) xmax = water[n][1];
if ( water[n][2] < ymin ) ymin = water[n][2];
if ( water[n][2] > ymax ) ymax = water[n][2];
n++;
}
n = 1;
while ( water[n][0] == 1 || water[n][0] == 2 ) {
p2.x = water[n][1];
p2.y = water[n][2];
p2.z = 0;
// System.out.println( "p1 " + p1.x + ", " + p1.y + ", " + p1.z );
// System.out.println( "p2 " + p2.x + ", " + p2.y + ", " + p2.z );
if ( water[n][0] == water[n-1][0] ) {
v = findDistanceFromLine( p1, p2, p );
// System.out.println( "v " + v.x + ", " + v.y + ", " + v.z + " : " + v.vx );
if ( inside( v.vx, v.x ) ) {
if ( v.y > 0 ) {
if ( nphaseline == 1 ) {
System.out.println( " solid/liquid " );
state = 10;
} else {
System.out.println( " solid " );
state = 0;
break;
}
nphaseline++;
} else {
if ( nphaseline == 1 ) {
System.out.println( " vapour " );
state = 2;
break;
} else {
System.out.println( " liquid " );
state = 1;
break;
}
}
}
} else {
System.out.println( "skip " );
}
p1.x = p2.x;
p1.y = p2.y;
p1.z = p2.z;
n++;
// System.out.println( "line " + n );
}
System.out.println( " state = " + state );
dataValid( tK, xmax, xmin );
dataValid( pPa, ymax, ymin );
return( state );
}
// find distance from line b0 - b1 to b2
public static Point3d findDistanceFromLine( Point3d b0, Point3d b1, Point3d b2 ) {
double angle1, angle2, angle3, distance = 0;
Point3d v0 = new Point3d();
Point3d v1 = new Point3d();
// transform 3 tracked locations into 2 vectors on x-y plane
// First transform around z axis
v0.x = b1.x - b0.x;
v0.y = b1.y - b0.y;
v0.z = b1.z - b0.z;
angle1 = Math.atan2(v0.y, v0.x);
v0 = transformAroundZaxis( v0, -angle1 );
v1.x = b2.x - b0.x;
v1.y = b2.y - b0.y;
v1.z = b2.z - b0.z;
v1 = transformAroundZaxis( v1, -angle1 );
// Then transform around y axis
angle2 = Math.atan2(v0.z, v0.x);
v0 = transformAroundYaxis( v0, -angle2 );
v1 = transformAroundYaxis( v1, -angle2 );
// And finally transform around x axis
// angle3 = Math.atan2(v1.z, v1.y);
// v1 = transformAroundXaxis( v1, -angle3 );
// distance = v1.y;
// System.out.println( " " + v0.x + " + " + v0.y + " + " + v0.z );
// v0 lies along x axis, and so v0.x is the length of the line element vector
// and if v1 is > 0 oe < v0.x, it is adjacent to v0, and either above it (v1.y > 0)
// or below it (v1.y < 0)
v1.vx = v0.x;
return( v1 );
}
// return true if x1 <= x0 and x1 >= 0
public static boolean inside( double x0, double x1 ) {
boolean inside = false;
if ( x1 >= 0 ) {
if ( x1 <= x0 ) {
inside = true;
}
}
return( inside );
}
public static boolean dataValid( double v, double vmax, double vmin ) {
boolean valid = false;
if ( v <= vmax && v >= vmin ) {
valid = true;
} else {
System.out.println( "Value " + v + " is invalid, outside range " + vmin + "--" + vmax);
}
return( valid );
}
// transform a point around the x-axis through xangle (june 2013)
public static Point3d transformAroundXaxis( Point3d m, double xangle ) {
Point3d mt = new Point3d();
mt.x = m.x;
mt.y = m.y * Math.cos(xangle) - m.z * Math.sin(xangle);
mt.z = m.z * Math.cos(xangle) + m.y * Math.sin(xangle);
mt.vx = m.vx;
mt.vy = m.vy * Math.cos(xangle) - m.vz * Math.sin(xangle);
mt.vz = m.vz * Math.cos(xangle) + m.vy * Math.sin(xangle);
return ( mt );
}
// transform a point around the y-axis through xangle (june 2013)
public static Point3d transformAroundYaxis( Point3d m, double xangle ) {
Point3d mt = new Point3d();
mt.y = m.y;
mt.x = m.x * Math.cos(xangle) - m.z * Math.sin(xangle);
mt.z = m.z * Math.cos(xangle) + m.x * Math.sin(xangle);
mt.vy = m.vy;
mt.vx = m.vx * Math.cos(xangle) - m.vz * Math.sin(xangle);
mt.vz = m.vz * Math.cos(xangle) + m.vx * Math.sin(xangle);
return ( mt );
}
// transform a point around the z-axis through xangle (june 2013)
public static Point3d transformAroundZaxis( Point3d m, double xangle ) {
Point3d mt = new Point3d();
mt.z = m.z;
mt.y = m.y * Math.cos(xangle) + m.x * Math.sin(xangle);
mt.x = m.x * Math.cos(xangle) - m.y * Math.sin(xangle);
mt.vz = m.vz;
mt.vy = m.vy * Math.cos(xangle) + m.vx * Math.sin(xangle);
mt.vx = m.vx * Math.cos(xangle) - m.vy * Math.sin(xangle);
return ( mt );
}
}
class Point3d {
double x;
double y;
double z;
double vx;
double vy;
double vz;
double ax;
double ay;
double az;
public Point3d() {
}
public Point3d( double x, double y, double z ) {
this.x = x;
this.y = y;
this.z = z;
}
}
// class Average calculates averages of values
class Average{
double runningTotal;
int nValues;
int averageOver;
double currentAverage;
public Average(){
}
public Average( int avOver ) {
this.runningTotal = 0;
this.nValues = 0;
this.averageOver = avOver;
}
// returns true if a new average value is available
boolean addNewValue( double value ) {
boolean newAverageAvailable = false;
this.runningTotal += value;
this.nValues++;
if ( nValues == averageOver ) {
currentAverage = this.runningTotal / (double)this.averageOver;
this.runningTotal = 0;
this.nValues = 0;
newAverageAvailable = true;
}
return( newAverageAvailable );
}
}
class PseudoMap{
// default world map is of present day Earth with all continents in N hemisphere
// world map of fractions of eqch latitude that are ocean. land, and snow
// and fraction of sphere surface area covered by each latitude
double pseudoMap[][] = {
// ocean land snow fAreaSphere Lat
{ 0.1, 0.4, 0.8, }, // albedos
{ 0.4, 0.6, 0, 0.043577872, }, // 2.5 40% ocean, 60% land
{ 0.4, 0.6, 0, 0.043246217, }, // 7.5
{ 0.4, 0.6, 0, 0.042585433, }, // 12.5
{ 0.4, 0.6, 0, 0.041600548, }, // 17.5
{ 0.4, 0.6, 0, 0.040299058, }, // 22.5
{ 0.4, 0.6, 0, 0.03869087, }, // 27.5
{ 0.4, 0.6, 0, 0.036788218, }, // 32.5
{ 0.4, 0.6, 0, 0.034605585, }, // 37.5
{ 0.4, 0.6, 0, 0.032159586, }, // 42.5
{ 0.4, 0.6, 0, 0.02946883, }, // 47.5
{ 0.4, 0.6, 0, 0.0265538, }, // 52.5
{ 0.4, 0.6, 0, 0.02343668, }, // 57.5
{ 0.4, 0.6, 0, 0.020141192, }, // 62.5
{ 0.4, 0.6, 0, 0.016692417, }, // 67.5
{ 0.4, 0.6, 0, 0.013116603, }, // 72.5
{ 0.4, 0.6, 0, 0.009440963, }, // 77.5
{ 0.4, 0.6, 0, 0.0056934725, }, // 82.5
{ 0, 0, 1.0, 0.0019026509, }, // 87.5 snow at north pole
{ 1.0, 0, 0, 0.043577872, }, // -2.5 100% ocean
{ 1.0, 0, 0, 0.043246217, }, // -7.5
{ 1.0, 0, 0, 0.042585433, }, //-12.5
{ 1.0, 0, 0, 0.041600548, }, //-17.5
{ 1.0, 0, 0, 0.040299058, }, //-22.5
{ 1.0, 0, 0, 0.03869087, }, //-27.5
{ 1.0, 0, 0, 0.036788218, }, //-32.5
{ 1.0, 0, 0, 0.034605585, }, //-37.5
{ 1.0, 0, 0, 0.032159586, }, //-42.5
{ 1.0, 0, 0, 0.02946883, }, //-47.5
{ 1.0, 0, 0, 0.0265538, }, //-52.5
{ 1.0, 0, 0, 0.02343668, }, //-57.5
{ 1.0, 0, 0, 0.020141192, }, //-62.5
{ 1.0, 0, 0, 0.016692417, }, //-67.5
{ 1.0, 0, 0, 0.013116603, }, //-72.5
{ 1.0, 0, 0, 0.009440963, }, //-77.5
{ 1.0, 0, 0, 0.0056934725, }, //-82.5
{ 0, 0, 1.0, 0.0019026509, } //-87.5 snow at south pole
};
double oceanAlbedo;
double landAlbedo;
double snowAlbedo;
public PseudoMap() {
oceanAlbedo = pseudoMap[0][0];
landAlbedo = pseudoMap[0][1];
snowAlbedo = pseudoMap[0][2];
System.out.println( "pseudoMap albedos " + oceanAlbedo + ", " + landAlbedo + ", " + snowAlbedo );
}
// aLat is surface area of latitude circle, a Sphere is surface area of sphere
double setLatitudeAreaFraction( int i, double aLat, double aSphere ) {
pseudoMap[i+1][3] = aLat / aSphere;
return( pseudoMap[i+1][3] );
}
// glaciate planet from north pole down to limitLatitude
void glaciatePseudoMap( double nLimitLatitude ) {
for ( int y=0; y<18; y++ ) {
if ( y > nLimitLatitude ) {
pseudoMap[y+1][0] = 0; // no ocean
pseudoMap[y+1][1] = 0; // no land
pseudoMap[y+1][2] = 1; // all snow
}
}
}
// glaciate planet from north pole down to limitLatitude
void glaciateLandPseudoMap( int nLimitLatitude ) {
for ( int y=0; y<18; y++ ) {
if ( y > nLimitLatitude ) {
pseudoMap[y+1][1] = 0; // no land
pseudoMap[y+1][2] = 1.0 - pseudoMap[y+1][0]; // all snow
}
}
}
void glaciateLand( int nLat ) {
pseudoMap[nLat+1][1] = 0; // no land
pseudoMap[nLat+1][2] = 1.0 - pseudoMap[nLat+1][0]; // all snow
}
void deglaciateLand( int nLat ) {
pseudoMap[nLat+1][1] = 1.0 - pseudoMap[nLat+1][0]; // all land
pseudoMap[nLat+1][2] = 0.0; // no snow
}
void printPseudoMap() {
for ( int y=0; y<18; y++ ) {
System.out.println( "nLatitude " + y + ": " + pseudoMap[y+1][0] + ", " + pseudoMap[y+1][1] + ", " + pseudoMap[y+1][2] + ", " + pseudoMap[y+1][3] + ", " );
}
}
// simple 3-colour rectangular world map of latitude lines
// paints from latitude +90 to -90
// (completely untested)
public void paint(Graphics g ) {
int x0 = 0;
int x1, x2, x3 = 72;
int latIndex;
g.setPaintMode();
for ( int y=0; y<36; y++ ) {
latIndex = y;
if ( y < 18 ) latIndex = 17 - y;
x1 = (int)( (double)x3 * pseudoMap[latIndex][1] ); // land
x2 = (int)( (double)x3 * pseudoMap[latIndex][2] ); // snow
g.setColor( Color.orange ); // land
g.drawLine(x0, y, x1, y);
g.setColor( Color.white ); // snow
g.drawLine(x1, y, x2, y);
g.setColor( Color.cyan ); // ocean
g.drawLine(x2, y, x3, y);
}
}
}
// daily values, mostly of air temperatures
class DailyValues {
double dailyValue[][] = new double[2][366];
double meanDailyValue;
double vsum = 0;
double albedo;
public DailyValues() {
dailyValue[1][0] = 0;
}
void setDailyValue( double v, int i ) {
dailyValue[0][i] = v;
vsum+= v;
dailyValue[1][i] = vsum;
}
double getDailyValue( int i ) {
return( dailyValue[0][i] );
}
// set daily value and calculate mean daily value
void setMeanDailyValue( double v, int i ) {
dailyValue[0][i] = v;
// assumes this method is repeatedly being called with i from 0 to 364
if ( i == 0 ) {
meanDailyValue = v;
} else if ( i == 364 ) {
meanDailyValue += v;
meanDailyValue = meanDailyValue / 365.0;
} else {
meanDailyValue += v;
}
}
// calculate mean daily value
double getMeanDailyValue() {
meanDailyValue = 0;
for ( int i=1; i<=365; i++ ) {
meanDailyValue += dailyValue[0][i];
}
meanDailyValue = meanDailyValue / 365.0;
return( meanDailyValue );
}
// get mean value between startDay and endDay
double getMeanDailyValue( int startDay, int endDay ) {
// System.out.println( "startDay" + startDay + " endDay " + endDay );
double period = (double)( endDay - startDay );
double vMean = ( this.dailyValue[1][endDay] - this.dailyValue[1][startDay] ) / period;
return( vMean );
}
void setDailySouthernHemisphereValue( double v, int i ) {
// southern hemisphere is offset 365/2 days
int is = i + 182;
if ( is >= 365 ) is = is - 365;
dailyValue[0][is] = v;
}
void setMeanDailySouthernHemisphereValue( double v, int i ) {
// southern hemisphere is offset 365/2 days
int is = i + 182;
if ( is >= 365 ) is = is - 365;
dailyValue[0][is] = v;
// assumes this method is repeatedly being called with i from 1 to 365
if ( is == 1 ) {
meanDailyValue = v;
} else if ( is == 365 ) {
meanDailyValue += v;
meanDailyValue = meanDailyValue / 365.0;
} else {
meanDailyValue += v;
}
}
void printDailyValues( String name) {
System.out.println( name );
for ( int i=1; i<=365; i++ ) {
System.out.print( (float)dailyValue[0][i] + " " );
if ( i % 5 == 0 ) System.out.println();
}
System.out.println();
}
}
subglacial 