Main Code

Static Forces (E & G) Between Particles in Debye Cells
Main Code
<?php
											 
											/*
												Main Code
												© Charles Chandler
												http://qdl.scs-inc.us/?top=11048
											*/
											 
											/*
												This calculates the gravity and electric forces between the core of a Debye cell and all of the
												particles in 56 other neighboring cells (cores and sheaths). Since gravity is purely attractive,
												we know that there will be a net body force pulling all of the cells together. The question is:
												will there be any net electric force, and if so, how much will it be, and will it be attractive
												or repulsive?
												 
												So we used another program to generate a random set of points representing the distribution of
												particles in a Debye sheath. We read this into the $ions array.
												 
												http://127.0.0.1/SCS/Other/QuickDisclosure/2ndParty/Images/Charles/Plasmas/DebyeForces.php
											*/
											 
											// Options
											define('offsetSheath',   0); // whether or not to offset the Debye sheaths from the dust grains
											define('doScatterPlots', 0); // plot the nuclei and +ions, to visually verify the spacing
											define('writeTextFile',  1); // write the x/y to a tab-delimited file
											 
											// Benchmarking
											@set_time_limit(0);
											$begTime = microtime(true);
											// elapsed: 58.87 for 5 iterations
											// elapsed: 281.23 for 20 iterations
											 
											// Includes
											include('GraphPoints.php');
											 
											// Constants
											define('e', 1.60217657 * pow(10, -19)); // elementary charge, in coulombs
											define('p', 1.67262178 * pow(10, -27)); // mass of proton, in kilograms
											define('neutralsPerIon', pow(10,  14)); // degree of ionization
											define('ElectricConstant', 8.98755 * pow(10,   9));
											define('GravityConstant',  6.67408 * pow(10, -11));
											$nucMass = 150 * neutralsPerIon * p;
											$ionMass =       neutralsPerIon * p;
											$scatterPlots = NULL;
											$newtonsPerKg = array();
											 
											// Functions
											function Pt($x, $y, $z = 0) {
											// just a shorthand for building the array
											return array('x' => $x, 'y' => $y, 'z' => $z);
											}
											 
											function PtOp($pt, $operator, $mixed) {
											if ($operator == '*') {
											if (is_numeric($mixed)) {
											return Pt($pt['x'] * $mixed, $pt['y'] * $mixed, $pt['z'] * $mixed);
											}
											elseif (is_array($mixed)) {
											return Pt($pt['x'] * $mixed['x'], $pt['y'] * $mixed['y'], $pt['z'] * $mixed['z']);
											}
											}
											elseif ($operator == '/') {
											if (is_numeric($mixed)) {
											return Pt($pt['x'] / $mixed, $pt['y'] / $mixed, $pt['z'] / $mixed);
											}
											elseif (is_array($mixed)) {
											return Pt($pt['x'] / $mixed['x'], $pt['y'] / $mixed['y'], $pt['z'] / $mixed['z']);
											}
											}
											elseif ($operator == '+') {
											if (is_numeric($mixed)) {
											return Pt($pt['x'] + $mixed, $pt['y'] + $mixed, $pt['z'] + $mixed);
											}
											elseif (is_array($mixed)) {
											return Pt($pt['x'] + $mixed['x'], $pt['y'] + $mixed['y'], $pt['z'] + $mixed['z']);
											}
											}
											elseif ($operator == '-') {
											if (is_numeric($mixed)) {
											return Pt($pt['x'] - $mixed, $pt['y'] - $mixed, $pt['z'] - $mixed);
											}
											elseif (is_array($mixed)) {
											return Pt($pt['x'] - $mixed['x'], $pt['y'] - $mixed['y'], $pt['z'] - $mixed['z']);
											}
											}
											}
											 
											function Gravity_n($m1_kg, $m2_kg, $dist_m) {
											// Given two masses (in kilograms), and the distance between them
											// (in meters), this returns the gravitational force (in Newtons).
											return (GravityConstant * $m1_kg * $m2_kg) / ($dist_m * $dist_m);
											}
											 
											function Electricity_n($num_e_1, $num_e_2, $dist_m) {
											// Given two quantities of elementary charges, and the distance between
											// them (in meters), this returns the electric force (in Newtons).
											$c1 = $num_e_1 * e; // convert elementary charges to coulombs
											$c2 = $num_e_2 * e; // convert elementary charges to coulombs
											return (ElectricConstant * $c1 * $c2) / ($dist_m * $dist_m);
											}
											 
											function ScientificNotation($num, $prec) {
											$minus = ($num < 0) ? '&minus;' : '';
											$num = abs($num);
											$exp = NULL;
											if ($num) {
											$exp = round(log($num, 10));
											if (abs($exp) > 2) {
											$num = $num / pow(10, $exp);
											if ($num < 1) { $num *= 10; $exp -= 1; }
											if ($exp < 0) $exp = '&minus;'.(abs($exp));
											$exp = '&thinsp;&times;&thinsp;10<sup>'.$exp.'</sup>';
											} else $exp = NULL;
											}
											return $minus.round($num, $prec).$exp;
											}
											 
											function Distance3D($pt1, $pt2) {
											$dX = ($pt2['x'] - $pt1['x']);
											$dY = ($pt2['y'] - $pt1['y']);
											$dZ = ($pt2['z'] - $pt1['z']);
											return sqrt(($dX * $dX) + ($dY * $dY) + ($dZ * $dZ));
											}
											 
											// M A I N
											echo trim('
												<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
												<html>
												<head>
												<link rel="stylesheet" type="text/css" href="../../../../style.css">
												</head>
												<body>
											');
											 
											// I O N S   I N   S H E A T H
											$ions = array();
											foreach(explode("\n", file_get_contents('DebyeIons.txt')) as $line) {
											if ($line) {
											list($x, $y, $z) = explode(chr(9), $line);
											$ions[] = Pt($x, $y, $z);
											}
											}
											 
											// C L O S E S T   P A C K E D   A R R A N G E M E N T   O F   C E L L S
											$cellCenterPts = array();
											foreach(explode("\n", file_get_contents('DebyeHCP.txt')) as $line) {
											if ($line) {
											list($x, $y, $z) = explode(chr(9), $line);
											$cellCenterPts[] = Pt($x, $y, $z);
											}
											}
											 
											// C A L C U L A T E   G   A N D   E
											// Now calculate the gravitational and electric forces, between the origin and the particles in each
											// of the other Debye cells. We have the centers of each of the cells in the $cellCenterPts array, so
											// first we find the g & e forces between the negatively charged Debye nucleus at that location and
											// the test mass/charge (i.e., $origin). Then we can take our ion points, shift them so that they
											// surround the $cellCenterPts point of the cell being calculated, and then find the g & e forces from
											// the test mass/charge to each ion in the other cell.
											echo '
												<table class="standard">
												<tr>
												<th>cell<br />spacing<br />(m)</th>
												<th>N<sub>E</sub></th>
												<th>N<sub>G</sub></th>
												<th>e/g</th>
												</tr>
											';
											$graphPts = array();
											$origin = Pt(0, 0, 0); // the coords of the test charge
											$iterations = 0;
											// Recalculate for a range of distances between each Debye cell. To do this, we
											// just multiply the points in the $cellCenterPts array by a scalar before using them.
											// They start out being 2 meters apart in the data file, and we want to run
											// numbers up to 10 meters apart, so the upper limit is 5.
											for ($cenToCenScalar = 1; $cenToCenScalar <= 5.01; $cenToCenScalar += .1) {
											$iterations ++;
											$g = $e = 0;
											$plotPts = array(); // for plotting the points to be calculated, to visually verify
											foreach($cellCenterPts as $cenPt) {
											// nucleus to nucleus
											// Note that (0, 0, 0) is our test charge, and we
											// don't want to compare it to itself, so we check for it.
											$scaledCenPt = PtOp($cenPt, '*', $cenToCenScalar);
											if ($cenPt != $origin) {
											$dist = Distance3D($origin, $scaledCenPt);
											$g += Gravity_n($nucMass, $nucMass, $dist);
											$e -= Electricity_n(150, 150, $dist);
											if ($cenPt['z'] == 3.26) { // create a fat point for the nucleus
											$plotPts[] = Pt($scaledCenPt['x'] + 0.00, $scaledCenPt['y'] + 0.00);
											$plotPts[] = Pt($scaledCenPt['x'] + 0.00, $scaledCenPt['y'] + 0.01);
											$plotPts[] = Pt($scaledCenPt['x'] + 0.01, $scaledCenPt['y'] + 0.00);
											$plotPts[] = Pt($scaledCenPt['x'] + 0.01, $scaledCenPt['y'] + 0.01);
											}
											 
											// Throw in the forces between each of the ions attached to the central
											// cell and all of the other particles.
											foreach($ions as $originIonPt) {
											$dist = Distance3D($originIonPt, $scaledCenPt);
											$g += Gravity_n($ionMass, $nucMass, $dist);
											$e += Electricity_n(1, 150, $dist);
											 
											foreach($ions as $ionPt) {
											$tmpIonPt = PtOp($scaledCenPt, '+', $ionPt);
											$dist = Distance3D($originIonPt, $tmpIonPt);
											$g += Gravity_n($ionMass, $ionMass, $dist);
											$e -= Electricity_n(1, 1, $dist);
											}
											}
											}
											 
											// nucleus to +ions
											if (($cenPt != $origin) or (offsetSheath)) {
											// Since the cell-to-cell spacing is 2 m in the data file, and since we want to offset
											// the sheaths half the distance to the next cell, to put them right between the two
											// nuclei, that would be (($cenToCenScalar * 2) / 2), which is just the $cenToCenScalar.
											$sheathOffset = (offsetSheath) ? $cenToCenScalar : 0;
											foreach($ions as $ionPt) {
											$tmpIonPt = PtOp($scaledCenPt, '+', $ionPt);
											$tmpIonPt['x'] += $sheathOffset;
											if ($cenPt['z'] == 3.26) $plotPts[] = $tmpIonPt;
											$dist = Distance3D($origin, $tmpIonPt);
											$g += Gravity_n($nucMass, $ionMass, $dist);
											$e += Electricity_n(150, 1, $dist);
											}
											}
											}
											echo '
												<tr>
												<td>'.round($cenToCenScalar * 2, 2).'</td>
												<td align="right">'.ScientificNotation($e, 2).'</td>
												<td align="right">'.ScientificNotation($g, 2).'</td>
												<td align="right">'.round($e / $g, 2).'</td>
												</tr>
											';
											// Store the newtons per kilogram, so that we can
											// find the average force causing the collapse, at
											// least for 4/5 the way, so we can estimate the time.
											if ($cenToCenScalar > 1) $newtonsPerKg[] = (($e + $g) / $nucMass);
											 
											if (doScatterPlots) {
											$fileName = 'Spacing_'.$cenToCenScalar.'.png';
											$scatterPlots .= '<p>'.$fileName.'</p>';
											$scatterPlots .= '<p>'.GraphPoints($plotPts, 'Blank.png', $fileName, false).'</p>';
											}
											 
											// If the sheaths are still attached to the nuclei, we're getting repulsion,
											// not attraction, but we'd rather see all positive numbers in the graph.
											if (!offsetSheath) $e = -$e;
											$graphPts[] = Pt($cenToCenScalar * 2, $e / $g);
											}
											$endTime = microtime(true);
											$elapsed = $endTime - $begTime;
											 
											echo '</table>';
											echo '<br />';
											echo '<br />';
											echo '<p>'.GraphPoints($graphPts, 'Blank.png', 'Results.png', true).'</p>';
											echo '<br />';
											echo '<br />';
											echo '<p>'.(array_sum($newtonsPerKg) / count($newtonsPerKg)).'</p>';
											echo '<br />';
											echo '<br />';
											echo $scatterPlots;
											echo '<br />';
											echo 'elapsed: '.round($elapsed, 2).' for '.$iterations.' iterations';
											echo '<br />';
											echo '</body></html>';
											 
											if (writeTextFile) {
											$results = '';
											foreach($graphPts as $pt) $results .= ($pt['x']).chr(9).($pt['y']).chr(13);
											file_put_contents('Results.txt', $results);
											}
											 
										?>
← PREV	Powered by Quick Disclosure Lite © 2010~2021 SCS-INC.US	UP ↑