[GRASS-SVN] r56945 - grass-addons/grass7/raster/r.crater
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Jun 28 04:26:00 PDT 2013
Author: ychemin
Date: 2013-06-28 04:26:00 -0700 (Fri, 28 Jun 2013)
New Revision: 56945
Modified:
grass-addons/grass7/raster/r.crater/crater.c
grass-addons/grass7/raster/r.crater/crater.h
grass-addons/grass7/raster/r.crater/main.c
grass-addons/grass7/raster/r.crater/r.crater.html
Log:
basic functionality added
Modified: grass-addons/grass7/raster/r.crater/crater.c
===================================================================
--- grass-addons/grass7/raster/r.crater/crater.c 2013-06-28 07:04:08 UTC (rev 56944)
+++ grass-addons/grass7/raster/r.crater/crater.c 2013-06-28 11:26:00 UTC (rev 56945)
@@ -13,6 +13,8 @@
* solid_rock = boolean, is target made of solid rock?
*/
+#define PI 3.1415927
+
/* 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);
@@ -20,37 +22,88 @@
}
/* Kinetic Energy 1/2*m*v^2 (W in Melosh) */
-double kinetic_energy(double Mi,double Vi){
+double kinetic_energy(double r_proj, double L, double Vi){
+ double Mi = mass_impactor(r_proj,L);
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){
+/******************************************************************/
+/* Forward mode equations */
+/******************************************************************/
+/*Gault Scaling (Gault, 1974), equations 7.8.1a/b/c in Meloch*/
+double Gault_Dat(double W, double r_proj, double r_targ, double theta, int target_type){
double Dat;
- if (solid_rock){
+ if (target_type==3){ /*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);
+ 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);
+ }
}
} 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)
+ 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)*/
+/*Yield Scaling (Nordyke, 1962), equation 7.8.3 in Meloch*/
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)*/
+/*Pi-group Scaling (Schmidt and Holsapple, 1983), equation 7.8.4 in Meloch*/
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));
}
+/******************************************************************/
+/* Backward mode equations */
+/******************************************************************/
+/*Gault Scaling (Gault, 1974), equations 7.8.1a/b/c in Meloch*/
+double Gault_L(double Dat, double Vi, double r_proj, double r_targ, double theta, int target_type){
+ double L;
+ if (Dat >= 100){
+ /*for craters > 100m up to 1000 m + any kind of target material*/
+ L=pow(Dat/(0.27*pow(r_proj,1/6.0)*pow(r_targ,-1/2.0)*pow(1/3.0*PI*Vi*Vi,0.28)*pow(1/2.0,3.28)*pow(sin(theta),1/3.0)),1/3.28);
+ }
+
+ if (target_type==3){ /*Solid rock*/
+ if (Dat < 10){
+ /*for craters up to 10 m + solid target rock*/
+ L=pow(Dat/(0.015*pow(r_proj,1/6.0)*pow(r_targ,-1/2.0)*pow(1/3.0*PI*Vi*Vi,0.37)*pow(1/2.0,3.37)*pow(sin(theta),2/3.0)),1/3.37);
+ } else {
+ /*for craters up to 100 m + loose target rock or regolith*/
+ L=pow(Dat/(0.25*pow(r_proj,1/6.0)*pow(r_targ,-1/2.0)*pow(1/3.0*PI*Vi*Vi,0.29)*pow(1/2.0,3.29)*pow(sin(theta),1/3.0)),1/3.29);
+ }
+ } else {
+ if (Dat < 100){
+ /*for craters up to 100 m + loose target rock or regolith*/
+ L=pow(Dat/(0.25*pow(r_proj,1/6.0)*pow(r_targ,-1/2.0)*pow(1/3.0*PI*Vi*Vi,0.29)*pow(1/2.0,3.29)*pow(sin(theta),1/3.0)),1/3.29);
+ }
+ }
+ return(L);
+}
+/*Yield Scaling (Nordyke, 1962), equation 7.8.3 in Meloch*/
+double Yield_L(double Vi, double r_proj, double r_targ, double Dat){
+ return pow(Dat/(0.0133*pow(1/3.0*PI*Vi*Vi,1/3.4)*pow(1/2.0,3/3.4)+1.51*pow(r_proj,1/2.0)*pow(r_targ,-1/2.0)),3.4/6.4);
+}
+
+/*Pi-group Scaling (Schmidt and Holsapple, 1983), equation 7.8.4 in Meloch*/
+double Pi_L(double Vi, double r_proj, double r_targ, double Dat, double g){
+ return pow(Dat/(1.8*pow(r_proj,0.11)*pow(r_targ,-1/3.0)*pow(g,-0.22)*pow(1/3.0*PI*Vi*Vi,0.22)*pow(1/2.0,0.66)),1/0.79);
+}
+
+
+
+
+
+
Modified: grass-addons/grass7/raster/r.crater/crater.h
===================================================================
--- grass-addons/grass7/raster/r.crater/crater.h 2013-06-28 07:04:08 UTC (rev 56944)
+++ grass-addons/grass7/raster/r.crater/crater.h 2013-06-28 11:26:00 UTC (rev 56945)
@@ -11,11 +11,27 @@
* 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);
/* Kinetic Energy 1/2*m*v^2 (W in Melosh) */
-double kinetic_energy(double Mi,double Vi);
+double kinetic_energy(double r_proj, double L, double Vi);
+
+/******************************************************************/
+/* Forward mode equations */
+/******************************************************************/
/*Gault Scaling (Gault, 1974)*/
-double Gault_Dat(double W, double r_proj, double r_targ, double theta, bool solid_rock);
+double Gault_Dat(double W, double r_proj, double r_targ, double theta, int target_type);
/*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);
+
+/******************************************************************/
+/* Backward mode equations */
+/******************************************************************/
+/*Gault Scaling (Gault, 1974), equations 7.8.1a/b/c in Meloch*/
+double Gault_L(double Dat, double Vi, double r_proj, double r_targ, double theta, int target_type);
+/*Yield Scaling (Nordyke, 1962), equation 7.8.3 in Meloch*/
+double Yield_L(double Vi, double r_proj, double r_targ, double Dat);
+/*Pi-group Scaling (Schmidt and Holsapple, 1983), equation 7.8.4 in Meloch*/
+double Pi_L(double Vi, double r_proj, double r_targ, double Dat, double g);
Modified: grass-addons/grass7/raster/r.crater/main.c
===================================================================
--- grass-addons/grass7/raster/r.crater/main.c 2013-06-28 07:04:08 UTC (rev 56944)
+++ grass-addons/grass7/raster/r.crater/main.c 2013-06-28 11:26:00 UTC (rev 56945)
@@ -5,7 +5,6 @@
* AUTHOR(S): Yann Chemin - yann.chemin at gmail.com
* PURPOSE: Creates craters from meteorites
* or meteorites from craters
- * original code was in fortran77 from Meloch
*
* COPYRIGHT: (C) 2013 by the GRASS Development Team
*
@@ -88,13 +87,13 @@
input8 = G_define_standard_option(G_OPT_R_INPUT);
input8->key = "transient_crater_diameter";
- input8->description = _("Default mode: Name of transient crater diameter raster map [kg/m^3]");
+ input8->description = _("Default mode: Name of transient crater diameter raster map [m]");
input8->required = NO;
input9 = G_define_standard_option(G_OPT_R_INPUT);
input9->key = "final_crater_diameter";
- input9->description = _("Default mode: Name of final crater diameter raster map [kg/m^3]");
- input8->required = NO;
+ input9->description = _("Default mode: Name of final crater diameter raster map [m]");
+ input9->required = NO;
output = G_define_standard_option(G_OPT_R_OUTPUT);
output->description = _("Name for projectile size (default) or crater size (-c) or crater creation time (-t) raster map [m] or [s]");
@@ -134,7 +133,7 @@
result = output->answer;
/*Default Mode: Estimate projectile size from crater diameter*/
- int comptype = 0;
+ int mode = 0;
/*Check if return of duration of impact was requested*/
int return_time = 0;
@@ -152,17 +151,18 @@
/*Flagged Mode: Estimate crater diameter from projectile size*/
infd_L = Rast_open_old(idiameter, "");
inrast_L = Rast_allocate_d_buf();
- comptype = 1;/*Switch to pass to crater function for non default mode*/
+ mode = 1;/*Switch to pass to crater function for non default mode*/
/*Projectile Diameter Size*/
}else{
/*Default Mode: Estimate projectile size from crater diameter*/
- infd_Dt = Rast_open_old(tcrater_diameter_transient, "");
- inrast_Dt = Rast_allocate_d_buf();
- /*Transient Crater Diameter*/
-
- infd_Dfinal = Rast_open_old(tcrater_diameter_final, "");
- inrast_Dfinal = Rast_allocate_d_buf();
- /*If known, the final crater diameter*/
+ if(input8->answer){ /*Transient Crater Diameter*/
+ infd_Dt = Rast_open_old(tcrater_diameter_transient, "");
+ inrast_Dt = Rast_allocate_d_buf();
+ }
+ if(input9->answer){ /*If known, the final crater diameter*/
+ infd_Dfinal = Rast_open_old(tcrater_diameter_final, "");
+ inrast_Dfinal = Rast_allocate_d_buf();
+ }
}
/***************************************************/
@@ -211,6 +211,7 @@
DCELL d_L;
DCELL d_Dt;
DCELL d_Dfinal;
+ DCELL d_W;
G_percent(row, nrows, 2);
/* read input maps */
@@ -261,7 +262,26 @@
} else {
d_L = 0.0;
}
- 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);
+ if(mode==0){
+ /*Forward mode: known projectile details (L,r_proj,Vi)*/
+ /*Solid_rock or not (1 or 0) for Gault (flag2)*/
+ d_W = kinetic_energy(d_r_proj,d_L,d_v);
+ if(scaling_law==1) /*flag 2 is Gault scaling law*/
+ d_Dt = Gault_Dat(d_W,d_r_proj,d_r_targ,d_theta,d_ttype);
+ else if (scaling_law==2)/*flag3 is Yield scaling law*/
+ d_Dt = Yield_Dat(d_W,d_r_proj,d_r_targ,d_L);
+ else /*default is Pi-scaling*/
+ d_Dt = Pi_Dat(d_W,d_r_proj,d_r_targ,d_L);
+ } else {
+ /*Backward mode*/
+ if(scaling_law==1) /*flag 2 is Gault scaling law*/
+ d_L = Gault_L(d_Dt,d_v,d_r_proj,d_r_targ,d_theta,d_ttype);
+ else if (scaling_law==2)/*flag3 is Yield scaling law*/
+ d_L = Yield_L(d_v,d_r_proj,d_r_targ,d_Dt);
+ else /*default is Pi-scaling*/
+ d_L = Pi_L(d_v,d_r_proj,d_r_targ,d_Dt,d_g);
+ }
+ /*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, mode);*/
outrast[col] = d;
}
}
Modified: grass-addons/grass7/raster/r.crater/r.crater.html
===================================================================
--- grass-addons/grass7/raster/r.crater/r.crater.html 2013-06-28 07:04:08 UTC (rev 56944)
+++ grass-addons/grass7/raster/r.crater/r.crater.html 2013-06-28 11:26:00 UTC (rev 56945)
@@ -1,8 +1,25 @@
<h2>DESCRIPTION</h2>
<em>r.crater</em> This program estimates the size of a gravity dominated impact crater or the projectile that made it.
-Three different estimates are presented, but the pi-scaling method is currently considered the best!n
<p>
+<em>Forward mode</em> This mode needs to know the projectile details<br>
+L: projectile diameter (m)<br>
+r_proj: projectile density (kg/m^3)<br>
+Vi: Projectile velocity (km/s)<br>
+theta: projectile impact angle (degrees) for Gault scaling law (flag2)<br>
+Solid_rock or not (1 or 0) for Gault scaling law (flag2)<br>
+<p>
+<em>Backward mode</em>This mode needs to know the crater details<br>
+
+
+<h2>NOTES</h2>
+Gault scaling law saturates at craters 1000 Diameter Apparent Transient, and was essentially designed for regolith (Moon surface).
+
+<p>
+Below is explanation from the Meloch Fortran code (not included because of copyright)
+<p>
+Three different estimates are presented, but the pi-scaling method is currently considered the best!
+<p>
Impact conditions:
argv[1]: enter the impact velocity in km/sec
argv[2]: enter the impact angle in degrees
@@ -21,7 +38,7 @@
Mode 1, crater size
Mode 2, projectile size
<p>
-Mode 1: Estimate crater diameter from projectile size*/
+Mode 1: Estimate crater diameter from projectile size
Mode 1 case: Projectile descriptors:
argv[8]: enter the projectile diameter in m
<p>
More information about the grass-commit
mailing list