[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