Discussion of Code Used in Prediction Model
1.1 Predicting Carolina bay Orientations
1.1 Predicting Carolina bay Orientations
1.1.1 Bearing Adjustments Based on The Spherical Earth Surface
1.1.1 Bearing Adjustments Based on The Spherical Earth Surface
1.1.2 Bearing Adjustments Based on A Rotating Earth
1.1.2 Bearing Adjustments Based on A Rotating Earth
1.1.3 Bearing Adjustments Based On Latitude
1.1.3 Bearing Adjustments Based On Latitude
1.2 The Analytical Model
1.2 The Analytical Model
Our Carolina bay catalogue database spreadsheet was enhanced to generate predictions based on the spherical-earth, loft-time and latitude adjustments, generating a solution for each catalogued field’s predicted bearings while variables such as ejecta average and terminal velocity are adjusted. Charts are generated to correlate the predicted values with measured values, allowing heuristic tuning of the variables to identify best-fit solutions simultaneously across all bay fields simultaneously. The spreadsheet also generates catalogue-wide sets of GE KML to visualize in the GIS tool. The model is heuristically focused on the latitude and longitude of the proposed Saginaw crater’s three control points (NE, Centroid and SW).
1.2.1 Variables
Variable:
- The Latitude and longitude of the bay location being evaluated.
Fixed:
- Crater latitude and longitude parameters (selected as
Centroid: 43.68N, 84.82W,
Rampart SW: 43.724N, 84.944W,
Rampart NE: 44.624N, 82.659W)
- A reasonable average ground-vector velocity for the ejecta as heuristically derived from our analysis. Satisfactory simultaneous solutions range from 1km/sec to 5 km/sec, with 3km/sec selected as default based on sensitivity testing)
- Three variables to generate a reasonable terminal velocity V (currently 381 m/sec) using
V = sqrt ( (2 * W) / (Cd * r * A),
where
W is droplet weight of an 1-meter diameter droplet, (in Newtons)
a density variable r (selected as 2,000 kg/meter^3) ,
and A is the frontal area of that droplet.
Cd is the Coefficient of drag of the ejecta droplet (selected as 0.3).
- - An incidence angle for ground impact to resolve ground-plane velocities (selected as 45º from vertical). The resulting ground vector of terminal velocity is 269 m/sec with the default variables.
1.2.2 Loft Time Calculation and Earth Rotation Longitude Offset Value
Longitude offset in degrees = Loft Time (minutes) /4
Where Loft Time = Loft Distance / Average Ground Velocity (default 3 km/sec)
The Loft Distance between the crater centroid location and that of the bay is derived using the spherical law of cosines, applied in the following Java routine. All variables in radians:
private double GreatCircleDistance(double lat1, double lon1, double lat2, double lon2) {
/* resolve distance between forePoint and farPoint input in radian
* @return The distance in Kilometres */
double dLat = (lat2 - lat1);
double dLon = (lon2 - lon1);
double a = Math.sin(dLat / 2) * Math.sin(dLat / 2) + Math.cos(lat1)
* Math.cos(lat2) * Math.sin(dLon / 2) * Math.sin(dLon / 2);
return (earthRadius * 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a)));}
1.2.3 Determining Baseline (Static Earth) Bearing
private double GreatCircleBearing(double lat1, double lon1, double lat2,
double lon2) {
double dLon = (lon2 - lon1);
double y = Math.sin(dLon) * Math.cos(lat2);
double x = Math.cos(lat1) * Math.sin(lat2) - Math.sin(lat1)
* Math.cos(lat2) * Math.cos(dLon);
double Bearing = 180 + (Math.atan2(y, x) * convrd);
return (Bearing * convdr); }
1.2.4 Applying Latitude-difference Adjustments to Base Bearing
1.2.4 Applying Latitude-difference Adjustments to Base Bearing
Figure 2
Figure 2 depicts the trigonometry used to adjust the arrival bearing prediction based on the latitude-driven W➔E earth ground and atmosphere rotational velocities. In our specific proposal, a relevant set of W➔E velocities would be the Saginaw crater – rotating at 1,270 km/hr – and a generic ejecta field such as Bishopville – rotating 112 km/hr faster. This ground velocity differential is combined with the terminal velocity as demonstrated as in figure 2.
We first decompose the ejecta’s terminal velocity vector (termv) to yield its W➔E component, using the base bearing ( “_C” references the centroid calculation, also done for _SW and _NE):
EWcomponent_C = termv * Math.cos(baseBearing_C) + c_deltaV;
Next, the bearing is re-composed using ACOS and the adjusted W➔E component:
componentBearing_C = Math.acos( EWcomponent_C / termv);
Finally, the bearing is cast into the proper compass quadrant (lost in the ACOS trigonometry) using the previously extracted compass quadrant of the non-rotating Earth base arrival bearing:
PredictedBearing_C = ((componentBearing_C * convrd) + (latQuadrant_C * 90));
The full Java code set includes routines to deal with compass cardinal point crossings and modulo 360º adjustments required for conversions between trigonometric notation (0º ➔ 360º) and compass bearing spherical notation (-180º ➔ +180º).
The model uses the predicted bearings to generate Google Earth kml visualization elements (paths, overlays and placemarks) for display in the GIS viewer. A sample Java routine, which generates the predicted bearing as an ”arrow” overlay, is provided here:
public String bearingArrowKML(double lat , Double lon , Double rotationValue ){
/** Returns a sting of kml to place a Bearing Arrow at a bay site given lat, lon and rotation inputs in radians **/
double distance = 2 ; /* km for building latlon box */
double NELat = Math.asin(Math.sin(lat)
* Math.cos(distance / earthRadius) + Math.cos(lat)
* Math.sin(distance / earthRadius) * Math.cos(45 * convdr));
double NELon = lon + Math.atan2(Math.sin(45 * convdr)
* Math.sin(distance / earthRadius) * Math.cos(lat),
Math.cos(distance / earthRadius) - Math.sin(lat)
* Math.sin(NELat));
double SWLat = Math.asin(Math.sin(lat)
* Math.cos(distance / earthRadius) + Math.cos(lat)
* Math.sin(distance / earthRadius) * Math.cos(225 * convdr));
double SWLon = lon Math.atan2(Math.sin(225 * convdr)
* Math.sin(distance / earthRadius) * Math.cos(lat),
Math.cos(distance / earthRadius) - Math.sin(lat)
* Math.sin(SWLat));
Double rotationValueBox = (180 -rotationValue*convrd) % 360;
final String kmlA = "
final String kmlB = "
"Bearing Calculator V 3.0
© Cintos 2010 ]]>
"
"
"
final String kmlC = "
final String kmlD = "
final String kmlE = "
final String kmlF = "
final String kmlG = "
return (kmlA + elementName + kmlB + NELat*convrd + kmlC + SWLat * convrd + kmlD + NELon * convrd + kmlE + SWLon * convrd + kmlF + rotationValueBox + kmlG);}
An example “GroundOverlay” kml string generated by the above routine is shown below. The example will visualize the Bishopville bearing prediction:
- 1.3 Independent Testing of Hypothesis
Geological Research by Cintos is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
Based on a work at Cintos.org.
Permissions beyond the scope of this license may be available at http://cintos.org/about.html.