[GRASS-SVN] r56592 - grass-addons/grass7/raster/r.crater
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jun 4 09:01:32 PDT 2013
Author: ychemin
Date: 2013-06-04 09:01:32 -0700 (Tue, 04 Jun 2013)
New Revision: 56592
Added:
grass-addons/grass7/raster/r.crater/crater.h
Modified:
grass-addons/grass7/raster/r.crater/crater.c
grass-addons/grass7/raster/r.crater/main.c
Log:
Starting the recoding of crater.c, created header file for it
Modified: grass-addons/grass7/raster/r.crater/crater.c
===================================================================
--- grass-addons/grass7/raster/r.crater/crater.c 2013-06-04 15:36:35 UTC (rev 56591)
+++ grass-addons/grass7/raster/r.crater/crater.c 2013-06-04 16:01:32 UTC (rev 56592)
@@ -1,2 +1,56 @@
-/* Removed because of copyright */
-/* in reimplementation from the Chapter 7: Scaling of crater dimensions, from Melosh*/
+#include <stdio.h>
+#include <math.h>
+/*
+ * Public Domain
+ * Defining crater equations after chapter 7 of Melosh
+ * Dat = Diameter Apparent Transient
+ * r_proj = density projectile (impactor)
+ * r_targ = density target
+ * theta = impacting angle
+ * Vi = impacting velocity
+ * Mi = impactor mass
+ * L = impactor diameter
+ * solid_rock = boolean, is target made of solid rock?
+*/
+
+/* Convert density and diameter into mass of impactor (Mi) */
+double mass_impactor(double r_proj, double L){
+ double volume = (4/3.0)*PI*pow((L/2.0),3);
+ return(r_proj*volume);
+}
+
+/* Kinetic Energy 1/2*m*v^2 (W in Melosh) */
+double kinetic_energy(double Mi,double Vi){
+ return 0.5*Mi*Vi*Vi;
+}
+
+/*Gault Scaling (Gault, 1974)*/
+double Gault_Dat(double W, double r_proj, double r_targ, double theta, bool solid_rock){
+ double Dat;
+ if (solid_rock){
+ /*for craters up to 10 m + solid target rock*/
+ Dat=0.015*pow(r_proj,1/6.0)*pow(r_targ,-1/2.0)*pow(W,0.37)*pow(sin(theta),2/3.0);
+ if (Dat > 10){
+ /*for craters up to 100 m + loose target rock or regolith*/
+ Dat=0.25*pow(r_proj,1/6.0)*pow(r_targ,-1/2.0)*pow(W,0.29)*pow(sin(theta),1/3.0);
+ }
+ } else {
+ /*for craters up to 100 m + loose target rock or regolith*/
+ Dat=0.25*pow(r_proj,1/6.0)*pow(r_targ,-1/2.0)*pow(W,0.29)*pow(sin(theta),1/3.0)
+ if (Dat > 100){
+ /*for craters > 100m up to 1000 m + any kind of target material*/
+ Dat=0.27*pow(r_proj,1/6.0)*pow(r_targ,-1/2.0)*pow(W,0.28)*pow(sin(theta),1/3.0);
+ }
+ return(Dat);
+}
+
+/*Yield Scaling (Nordyke, 1962)*/
+double Yield_Dat(double W, double r_proj, double r_targ, double L){
+ return (0.0133*pow(W,1/3.4) + 1.51*pow(r_proj,1/2.0)*pow(r_targ,-1/2.0)*L);
+}
+
+/*Pi-group Scaling (Schmidt and Holsapple, 1983)*/
+double Pi_Dat(double W, double r_proj, double r_targ, double L, double g){
+ return (1.8*pow(r_proj,0.11)*pow(r_targ,-1/3.0)*pow(g,-0.22)*pow(L,0.13)*pow(W,0.22));
+}
+
Added: grass-addons/grass7/raster/r.crater/crater.h
===================================================================
--- grass-addons/grass7/raster/r.crater/crater.h (rev 0)
+++ grass-addons/grass7/raster/r.crater/crater.h 2013-06-04 16:01:32 UTC (rev 56592)
@@ -0,0 +1,21 @@
+/*
+ * Public Domain
+ * Defining crater equations after chapter 7 of Melosh
+ * Dat = Diameter Apparent Transient
+ * r_proj = density projectile (impactor)
+ * r_targ = density target
+ * theta = impacting angle
+ * Vi = impacting velocity
+ * Mi = impactor mass
+ * L = impactor diameter
+ * solid_rock = boolean, is target made of solid rock?
+*/
+
+/* Kinetic Energy 1/2*m*v^2 (W in Melosh) */
+double kinetic_energy(double Mi,double Vi);
+/*Gault Scaling (Gault, 1974)*/
+double Gault_Dat(double W, double r_proj, double r_targ, double theta, bool solid_rock);
+/*Yield Scaling (Nordyke, 1962)*/
+double Yield_Dat(double W, double r_proj, double r_targ, double L);
+/*Pi-group Scaling (Schmidt and Holsapple, 1983)*/
+double Pi_Dat(double W, double r_proj, double r_targ, double L, double g);
Modified: grass-addons/grass7/raster/r.crater/main.c
===================================================================
--- grass-addons/grass7/raster/r.crater/main.c 2013-06-04 15:36:35 UTC (rev 56591)
+++ grass-addons/grass7/raster/r.crater/main.c 2013-06-04 16:01:32 UTC (rev 56592)
@@ -33,15 +33,15 @@
struct History history; /*metadata */
char *result; /*output raster name */
- int infd_v, infd_theta, infd_rhotarget;
- int infd_g, infd_rhoproj, infd_ttyp, infd_ttype;
+ int infd_v, infd_theta, infd_r_targ;
+ int infd_g, infd_r_proj, infd_ttyp, infd_ttype;
int infd_L, infd_Dt, infd_Dfinal;/*Modes parameters*/
int outfd;
char *ivelocity, *iangle, *idensity, *idiameter; /*Impactor*/
char *tg, *ttype, *tdensity; /*Target*/
char *tcrater_diameter_transient, *tcrater_diameter_final; /*Target crater*/
- void *inrast_v, *inrast_theta, *inrast_rhotarget;
- void *inrast_g, *inrast_ttype, *inrast_rhoproj;
+ void *inrast_v, *inrast_theta, *inrast_r_targ;
+ void *inrast_g, *inrast_ttype, *inrast_r_proj;
void *inrast_L, *inrast_Dt, *inrast_Dfinal;
DCELL *outrast;
@@ -174,9 +174,9 @@
inrast_theta = Rast_allocate_d_buf();
/*theta = Impact angle in degrees*/
- infd_rhotarget = Rast_open_old(tdensity, "");
- inrast_rhotarget = Rast_allocate_d_buf();
- /*rhotarget = Target Density in kg/m^3*/
+ infd_r_targ = Rast_open_old(tdensity, "");
+ inrast_r_targ = Rast_allocate_d_buf();
+ /*r_targ = Target Density in kg/m^3*/
infd_g = Rast_open_old(tg, "");
inrast_g = Rast_allocate_d_buf();
@@ -186,9 +186,9 @@
inrast_ttype = Rast_allocate_d_buf();
/*targype = liqH2O=1 Loose_Sand=2 Competent_rock/saturated_soil=3*/
- infd_rhoproj = Rast_open_old(idensity, "");
- inrast_rhoproj = Rast_allocate_d_buf();
- /*rhoproj = Density of Projectile */
+ infd_r_proj = Rast_open_old(idensity, "");
+ inrast_r_proj = Rast_allocate_d_buf();
+ /*r_proj = Density of Projectile */
/***************************************************/
nrows = Rast_window_rows();
@@ -204,8 +204,8 @@
DCELL d;
DCELL d_v;
DCELL d_theta;
- DCELL d_rhotarget;
- DCELL d_rhoproj;
+ DCELL d_r_targ;
+ DCELL d_r_proj;
DCELL d_g;
DCELL d_ttype;
DCELL d_L;
@@ -216,8 +216,8 @@
/* read input maps */
Rast_get_d_row(infd_v, inrast_v, row);
Rast_get_d_row(infd_theta, inrast_theta, row);
- Rast_get_d_row(infd_rhotarget, inrast_rhotarget, row);
- Rast_get_d_row(infd_rhoproj, inrast_rhoproj, row);
+ Rast_get_d_row(infd_r_targ, inrast_r_targ, row);
+ Rast_get_d_row(infd_r_proj, inrast_r_proj, row);
Rast_get_d_row(infd_g, inrast_g, row);
Rast_get_d_row(infd_ttype, inrast_ttype, row);
if(flag->answer && input7->answer){
@@ -232,8 +232,8 @@
{
d_v = ((DCELL *) inrast_v)[col];
d_theta = ((DCELL *) inrast_theta)[col];
- d_rhotarget = ((DCELL *) inrast_rhotarget)[col];
- d_rhoproj = ((DCELL *) inrast_rhoproj)[col];
+ d_r_targ = ((DCELL *) inrast_r_targ)[col];
+ d_r_proj = ((DCELL *) inrast_r_proj)[col];
d_g = ((DCELL *) inrast_g)[col];
d_ttype = ((DCELL *) inrast_ttype)[col];
if(flag->answer){
@@ -244,12 +244,12 @@
}
if (Rast_is_d_null_value(&d_v) ||
Rast_is_d_null_value(&d_theta) ||
- Rast_is_d_null_value(&d_rhotarget) ||
- Rast_is_d_null_value(&d_rhoproj) ||
+ Rast_is_d_null_value(&d_r_targ) ||
+ Rast_is_d_null_value(&d_r_proj) ||
Rast_is_d_null_value(&d_g) ||
Rast_is_d_null_value(&d_ttype) ||
((d_ttype < 1) || (d_ttype > 3)) ||
- Rast_is_d_null_value(&d_rhoproj) ||
+ Rast_is_d_null_value(&d_r_proj) ||
((input7->answer)&&(Rast_is_d_null_value(&d_L))) ||
((input8->answer)&&(Rast_is_d_null_value(&d_Dt))) ||
((input9->answer)&&(Rast_is_d_null_value(&d_Dfinal))) )
@@ -261,7 +261,7 @@
} else {
d_L = 0.0;
}
- d = crater(d_v, d_theta, d_rhotarget, d_g, d_ttype, d_rhoproj, d_L, d_Dt, d_Dfinal, scaling_law, return_time, comptype);
+ d = crater(d_v, d_theta, d_r_targ, d_g, d_ttype, d_r_proj, d_L, d_Dt, d_Dfinal, scaling_law, return_time, comptype);
outrast[col] = d;
}
}
@@ -269,12 +269,12 @@
}
G_free(inrast_v);
G_free(inrast_theta);
- G_free(inrast_rhoproj);
+ G_free(inrast_r_proj);
G_free(inrast_g);
Rast_close(infd_v);
Rast_close(infd_theta);
- Rast_close(infd_rhotarget);
- Rast_close(infd_rhoproj);
+ Rast_close(infd_r_targ);
+ Rast_close(infd_r_proj);
Rast_close(infd_g);
G_free(outrast);
Rast_close(outfd);
More information about the grass-commit
mailing list