[GRASS-SVN] r31483 - in grass-addons/gipe: . i.evapo.SENAY
r.soiltex2prop
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu May 22 09:41:51 EDT 2008
Author: ychemin
Date: 2008-05-22 09:41:51 -0400 (Thu, 22 May 2008)
New Revision: 31483
Added:
grass-addons/gipe/r.soiltex2prop/
grass-addons/gipe/r.soiltex2prop/Makefile
grass-addons/gipe/r.soiltex2prop/description.html
grass-addons/gipe/r.soiltex2prop/main.c
grass-addons/gipe/r.soiltex2prop/prct2hf.c
grass-addons/gipe/r.soiltex2prop/prct2ksat.c
grass-addons/gipe/r.soiltex2prop/prct2porosity.c
grass-addons/gipe/r.soiltex2prop/vector_multiplication.c
Modified:
grass-addons/gipe/Makefile
grass-addons/gipe/i.evapo.SENAY/main.c
Log:
added new module converting soil texture to porosity (vol. fraction), saturated hydraulic conductivity (cm/h) and wetting front suction (cm), after Rawls et al., 1990.
Modified: grass-addons/gipe/Makefile
===================================================================
--- grass-addons/gipe/Makefile 2008-05-22 11:22:51 UTC (rev 31482)
+++ grass-addons/gipe/Makefile 2008-05-22 13:41:51 UTC (rev 31483)
@@ -85,6 +85,7 @@
r.rescale.eq \
r.series \
r.slope.aspect \
+ r.soiltex2prop \
r.statistics \
r.stats \
r.sum \
Modified: grass-addons/gipe/i.evapo.SENAY/main.c
===================================================================
--- grass-addons/gipe/i.evapo.SENAY/main.c 2008-05-22 11:22:51 UTC (rev 31482)
+++ grass-addons/gipe/i.evapo.SENAY/main.c 2008-05-22 13:41:51 UTC (rev 31483)
@@ -36,7 +36,7 @@
struct GModule *module;
struct Option *input1, *input2, *input3, *input4;
struct Option *input5, *input6, *input7, *input8, *input9;
- struct Option *input10, *input11, *input12;
+ struct Option *input10, *input11, *input12, *input13;
struct Option *output1, *output2;
struct Flag *flag2, *flag3;
@@ -51,10 +51,11 @@
int infd_lat, infd_doy, infd_tsw;
int infd_slope, infd_aspect;
int infd_tair, infd_e0;
+ int infd_ndvi;
int outfd1, outfd2;
char *albedo,*tempk,*dem,*lat,*doy,*tsw;
- char *slope,*aspect,*tair,*e0;
+ char *slope,*aspect,*tair,*e0, *ndvi;
double roh_w, e_atm;
int i=0,j=0;
@@ -63,6 +64,7 @@
void *inrast_doy, *inrast_tsw;
void *inrast_slope, *inrast_aspect;
void *inrast_tair, *inrast_e0;
+ void *inrast_ndvi;
DCELL *outrast1, *outrast2;
RASTER_MAP_TYPE data_type_output=DCELL_TYPE;
@@ -76,6 +78,7 @@
RASTER_MAP_TYPE data_type_aspect;
RASTER_MAP_TYPE data_type_tair;
RASTER_MAP_TYPE data_type_e0;
+ RASTER_MAP_TYPE data_type_ndvi;
/********************************/
/* Stats for Senay equation */
double t0dem_min=400.0,t0dem_max=200.0;
@@ -122,7 +125,6 @@
input7->key =_("roh_w");
input7->type = TYPE_DOUBLE;
input7->required = YES;
- input7->gisprompt =_("value, parameter");
input7->description=_("Value of the density of fresh water ~[1000-1020]");
input7->answer =_("1005.0");
@@ -142,7 +144,6 @@
input10->key =_("e_atm");
input10->type = TYPE_DOUBLE;
input10->required = NO;
- input10->gisprompt =_("value, parameter");
input10->description=_("Value of the apparent atmospheric emissivity (Bandara, 1998 used 0.845 for Sri Lanka)");
input10->guisection = _("Optional");
@@ -158,6 +159,11 @@
input12->description=_("Name of the Surface Emissivity map [-], use with -b");
input12->guisection = _("Optional");
+ input13 = G_define_standard_option(G_OPT_R_INPUT) ;
+ input13->key =_("ndvi");
+ input13->description=_("Name of the NDVI map [-]");
+ input13->answer =_("ndvi");
+
output1 = G_define_standard_option(G_OPT_R_OUTPUT) ;
output1->key =_("eta");
output1->description=_("Name of the output Actual ET layer");
@@ -189,9 +195,12 @@
roh_w = atof(input7->answer);
slope = input8->answer;
aspect = input9->answer;
- e_atm = atof(input10->answer);
+ if(input10->answer){
+ e_atm = atof(input10->answer);
+ }
tair = input11->answer;
e0 = input12->answer;
+ ndvi = input13->answer;
result1 = output1->answer;
result2 = output2->answer;
@@ -314,6 +323,17 @@
inrast_e0 = G_allocate_raster_buf(data_type_e0);
}
/***************************************************/
+ mapset = G_find_cell2 (ndvi, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("Cell file [%s] not found"), ndvi);
+ }
+ data_type_ndvi = G_raster_map_type(ndvi,mapset);
+ if ( (infd_ndvi = G_open_cell_old (ndvi,mapset)) < 0)
+ G_fatal_error(_("Cannot open cell file [%s]"), ndvi);
+ if (G_get_cellhd (ndvi, mapset, &cellhd) < 0)
+ G_fatal_error(_("Cannot read file header of [%s]"), ndvi);
+ inrast_ndvi = G_allocate_raster_buf(data_type_ndvi);
+ /***************************************************/
G_debug(3, "number of rows %d",cellhd.rows);
nrows = G_window_rows();
ncols = G_window_cols();
@@ -336,6 +356,7 @@
DCELL d_tempk;
DCELL d_dem;
DCELL d_t0dem;
+ DCELL d_ndvi;
G_percent(row,nrows,2);
if(G_get_raster_row(infd_albedo,inrast_albedo,row,data_type_albedo)<0)
G_fatal_error(_("Could not read from <%s>"),albedo);
@@ -343,6 +364,8 @@
G_fatal_error(_("Could not read from <%s>"),tempk);
if(G_get_raster_row(infd_dem,inrast_dem,row,data_type_dem)<0)
G_fatal_error(_("Could not read from <%s>"),dem);
+ if(G_get_raster_row(infd_ndvi,inrast_ndvi,row,data_type_ndvi)<0)
+ G_fatal_error(_("Could not read from <%s>"),ndvi);
/*process the data */
for (col=0; col < ncols; col++)
{
@@ -379,9 +402,21 @@
d_dem = (double) ((DCELL *) inrast_dem)[col];
break;
}
+ switch(data_type_ndvi){
+ case CELL_TYPE:
+ d_ndvi = (double) ((CELL *) inrast_ndvi)[col];
+ break;
+ case FCELL_TYPE:
+ d_ndvi = (double) ((FCELL *) inrast_ndvi)[col];
+ break;
+ case DCELL_TYPE:
+ d_ndvi = (double) ((DCELL *) inrast_ndvi)[col];
+ break;
+ }
if(G_is_d_null_value(&d_albedo)||
G_is_d_null_value(&d_tempk)||
- G_is_d_null_value(&d_dem)){
+ G_is_d_null_value(&d_dem)||
+ G_is_d_null_value(&d_ndvi)){
/* do nothing */
} else {
d_t0dem = d_tempk + 0.00649*d_dem;
@@ -389,7 +424,8 @@
/* do nothing */
} else {
//if(d_t0dem<t0dem_min){
- if(d_tempk<tempk_min&&d_albedo<0.1){
+ if(d_tempk<tempk_min&&
+ d_albedo<0.1){
t0dem_min=d_t0dem;
tempk_min=d_tempk;
//}else if(d_t0dem>t0dem_max){
@@ -419,6 +455,7 @@
DCELL d_aspect;
DCELL d_tair;
DCELL d_e0;
+ DCELL d_ndvi;
DCELL d_etpotd;
DCELL d_evapfr;
G_percent(row,nrows,2);
@@ -557,6 +594,7 @@
G_is_d_null_value(&d_lat)||
G_is_d_null_value(&d_doy)||
G_is_d_null_value(&d_tsw)||
+ G_is_d_null_value(&d_ndvi)||
((flag2->answer)&&G_is_d_null_value(&d_slope))||
((flag2->answer)&&G_is_d_null_value(&d_aspect))||
((flag3->answer)&&G_is_d_null_value(&d_tair))||
@@ -570,7 +608,15 @@
} else {
d_solar = solar_day(d_lat, d_doy, d_tsw );
}
- d_evapfr = evapfr_senay( tempk_max, tempk_min, d_tempk);
+ d_evapfr=evapfr_senay(tempk_max,tempk_min, d_tempk);
+ /*If water then no water stress*/
+ if(d_albedo<=0.1&&d_ndvi<=0.0){
+ d_evapfr=1.0;
+ }
+ /*some points are colder than selected low*/
+ if(d_evapfr>1.0){
+ d_evapfr=1.0;
+ }
if(result2){
outrast2[col] = d_evapfr;
}
@@ -580,7 +626,6 @@
d_rnetd = r_net_day(d_albedo,d_solar,d_tsw);
}
d_etpotd = et_pot_day(d_rnetd,d_tempk,roh_w);
- /*d_etpotd *= d_tsw;*/
d = d_etpotd * d_evapfr;
outrast1[col] = d;
}
Added: grass-addons/gipe/r.soiltex2prop/Makefile
===================================================================
--- grass-addons/gipe/r.soiltex2prop/Makefile (rev 0)
+++ grass-addons/gipe/r.soiltex2prop/Makefile 2008-05-22 13:41:51 UTC (rev 31483)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.soiltex2prop
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/gipe/r.soiltex2prop/description.html
===================================================================
--- grass-addons/gipe/r.soiltex2prop/description.html (rev 0)
+++ grass-addons/gipe/r.soiltex2prop/description.html 2008-05-22 13:41:51 UTC (rev 31483)
@@ -0,0 +1,24 @@
+<H2>DESCRIPTION</H2>
+
+<EM>r.soiltex2prop</EM> calculates the soil porosity (volume fraction), saturated hydraulic conductivity (cm/h) and wetting front suction (cm ; Green-Ampt) after Rawls et al, 1990.
+
+It takes input of soil texture proportions (sand, clay) in range [0.0-100.0].
+<H2>NOTES</H2>
+
+<H2>TODO</H2>
+
+
+<H2>SEE ALSO</H2>
+
+<em>
+<A HREF="r.uslek.html">r.uslek</A><br>
+</em>
+
+
+<H2>AUTHORS</H2>
+
+Yann Chemin, International Rice Research Institute, The Philippines<BR>
+
+
+<p>
+<i>Last changed: $Date: 2006/09/31 15:20:43 $</i>
Added: grass-addons/gipe/r.soiltex2prop/main.c
===================================================================
--- grass-addons/gipe/r.soiltex2prop/main.c (rev 0)
+++ grass-addons/gipe/r.soiltex2prop/main.c 2008-05-22 13:41:51 UTC (rev 31483)
@@ -0,0 +1,209 @@
+/****************************************************************************
+ *
+ * MODULE: r.soiltex2prop
+ * AUTHOR(S): Yann Chemin - ychemin at gmail.com
+ * PURPOSE: Transforms percentage of texture (sand/clay)
+ * into USDA 1951 (p209) soil texture classes and then
+ * into USLE soil erodibility factor (K) as an output
+ *
+ * COPYRIGHT: (C) 2008 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#define POLYGON_DIMENSION 50
+double point_in_triangle(double point_x, double point_y, double point_z, double t1_x, double t1_y, double t1_z, double t2_x, double t2_y, double t2_z, double t3_x, double t3_y, double t3_z);
+int prct2porosity(double sand_input, double clay_input);
+int prct2ksat(double sand_input, double clay_input);
+int prct2hf(double sand_input, double clay_input);
+
+int main(int argc, char *argv[])
+{
+ struct Cell_head cellhd; //region+header info
+ char *mapset; // mapset name
+ int nrows, ncols;
+ int row,col;
+
+ struct GModule *module;
+ struct Option *input1, *input2, *input3, *input4;
+ struct Option *output1, *output2, *output3;
+
+ struct Flag *flag1;
+ struct History history; //metadata
+
+ /************************************/
+ /* FMEO Declarations*****************/
+ char *name; // input raster name
+ char *result1; //output porosity raster name
+ char *result2; //output KSat raster name
+ char *result3; //output Hf raster name
+ //File Descriptors
+ int infd_psand, infd_pclay;
+ int outfd1, outfd2, outfd3;
+
+ char *psand,*pclay;
+
+ int i=0,j=0;
+
+ void *inrast_psand, *inrast_pclay;
+ DCELL *outrast1, *outrast2, *outrast3;
+ RASTER_MAP_TYPE data_type_output=DCELL_TYPE;
+ RASTER_MAP_TYPE data_type_psand;
+ RASTER_MAP_TYPE data_type_pclay;
+
+ /************************************/
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ module->keywords = _("soil, texture, porosity, Ksat, Hf");
+ module->description = _("Texture to soil properties (porosity, Ksat,Hf)");
+
+ /* Define the different options */
+ input1 = G_define_standard_option(G_OPT_R_INPUT) ;
+ input1->key = _("psand");
+ input1->description=_("Name of the Soil sand fraction map [0.0-1.0]");
+ input1->answer =_("psand");
+
+ input2 = G_define_standard_option(G_OPT_R_INPUT) ;
+ input2->key =_("pclay");
+ input2->description=_("Name of the Soil clay fraction map [0.0-1.0]");
+ input2->answer =_("pclay");
+
+ output1 = G_define_standard_option(G_OPT_R_OUTPUT) ;
+ output1->key =_("porosity");
+ output1->description=_("Name of the output porosity layer");
+ output1->answer =_("porosity");
+
+ output2 = G_define_standard_option(G_OPT_R_OUTPUT) ;
+ output2->key =_("ksat");
+ output2->description=_("Name of the output Ksat layer");
+ output2->answer =_("ksat");
+
+ output3 = G_define_standard_option(G_OPT_R_OUTPUT) ;
+ output3->key =_("hf");
+ output3->description=_("Name of the output hf layer");
+ output3->answer =_("hf");
+
+ /********************/
+ if (G_parser(argc, argv))
+ exit (EXIT_FAILURE);
+
+ psand = input1->answer;
+ pclay = input2->answer;
+
+ result1 = output1->answer;
+ result2 = output2->answer;
+ result3 = output3->answer;
+ /***************************************************/
+ mapset = G_find_cell2(psand, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("cell file [%s] not found"), psand);
+ }
+ data_type_psand = G_raster_map_type(psand,mapset);
+ if ( (infd_psand = G_open_cell_old (psand,mapset)) < 0)
+ G_fatal_error (_("Cannot open cell file [%s]"), psand);
+ if (G_get_cellhd (psand, mapset, &cellhd) < 0)
+ G_fatal_error (_("Cannot read file header of [%s])"), psand);
+ inrast_psand = G_allocate_raster_buf(data_type_psand);
+ /***************************************************/
+ mapset = G_find_cell2 (pclay, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("Cell file [%s] not found"), pclay);
+ }
+ data_type_pclay = G_raster_map_type(pclay,mapset);
+ if ( (infd_pclay = G_open_cell_old (pclay,mapset)) < 0)
+ G_fatal_error(_("Cannot open cell file [%s]"), pclay);
+ if (G_get_cellhd (pclay, mapset, &cellhd) < 0)
+ G_fatal_error(_("Cannot read file header of [%s]"), pclay);
+ inrast_pclay = G_allocate_raster_buf(data_type_pclay);
+ /***************************************************/
+ G_debug(3, "number of rows %d",cellhd.rows);
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ outrast1 = G_allocate_raster_buf(data_type_output);
+ outrast2 = G_allocate_raster_buf(data_type_output);
+ outrast3 = G_allocate_raster_buf(data_type_output);
+ /* Create New raster files */
+ if ( (outfd1 = G_open_raster_new(result1,data_type_output)) < 0)
+ G_fatal_error(_("Could not open <%s>"),result1);
+ if ( (outfd2 = G_open_raster_new(result2,data_type_output)) < 0)
+ G_fatal_error(_("Could not open <%s>"),result2);
+ if ( (outfd3 = G_open_raster_new(result3,data_type_output)) < 0)
+ G_fatal_error(_("Could not open <%s>"),result3);
+ /* Process pixels */
+ for (row = 0; row < nrows; row++)
+ {
+ DCELL d;
+ DCELL d_sand;
+ DCELL d_clay;
+ double tex[4] = {0.0,0.0,0.0,0.0};
+ G_percent(row,nrows,2);
+// printf("row = %i/%i\n",row,nrows);
+ /* read soil input maps */
+ if(G_get_raster_row(infd_psand,inrast_psand,row,data_type_psand)<0)
+ G_fatal_error(_("Could not read from <%s>"),psand);
+ if(G_get_raster_row(infd_pclay,inrast_pclay,row,data_type_pclay)<0)
+ G_fatal_error(_("Could not read from <%s>"),pclay);
+ /*process the data */
+ for (col=0; col < ncols; col++)
+ {
+ d_sand = ((DCELL *) inrast_psand)[col];
+ d_clay = ((DCELL *) inrast_pclay)[col];
+ if(G_is_d_null_value(&d_sand)||
+ G_is_d_null_value(&d_clay)){
+ G_set_d_null_value(&outrast1[col],1);
+ G_set_d_null_value(&outrast2[col],1);
+ G_set_d_null_value(&outrast3[col],1);
+ } else {
+ /************************************/
+ /* convert to porosity */
+ d = prct2porosity(d_sand, d_clay);
+ outrast1[col] = d;
+ d = prct2ksat(d_sand, d_clay);
+ outrast2[col] = d;
+ d = prct2hf(d_sand, d_clay);
+ outrast3[col] = d;
+ }
+ }
+ if (G_put_raster_row (outfd1, outrast1, data_type_output)<0)
+ G_fatal_error(_("Cannot write output raster file1"));
+ if (G_put_raster_row (outfd2, outrast2, data_type_output)<0)
+ G_fatal_error(_("Cannot write output raster file2"));
+ if (G_put_raster_row (outfd3, outrast3, data_type_output)<0)
+ G_fatal_error(_("Cannot write output raster file3"));
+ }
+
+ G_free (inrast_psand);
+ G_free (inrast_pclay);
+ G_close_cell (infd_psand);
+ G_close_cell (infd_pclay);
+
+ G_free (outrast1);
+ G_close_cell (outfd1);
+ G_free (outrast2);
+ G_close_cell (outfd2);
+ G_free (outrast3);
+ G_close_cell (outfd3);
+
+ G_short_history(result1, "raster", &history);
+ G_command_history(&history);
+ G_write_history(result1,&history);
+ G_short_history(result2, "raster", &history);
+ G_command_history(&history);
+ G_write_history(result2,&history);
+ G_short_history(result3, "raster", &history);
+ G_command_history(&history);
+ G_write_history(result3,&history);
+
+ exit(EXIT_SUCCESS);
+}
+
Added: grass-addons/gipe/r.soiltex2prop/prct2hf.c
===================================================================
--- grass-addons/gipe/r.soiltex2prop/prct2hf.c (rev 0)
+++ grass-addons/gipe/r.soiltex2prop/prct2hf.c 2008-05-22 13:41:51 UTC (rev 31483)
@@ -0,0 +1,716 @@
+#include<stdio.h>
+
+#define POLYGON_DIMENSION 50
+
+struct vector{
+ double sand;
+ double clay;
+ double silt;
+};
+
+
+// HF
+
+int prct2hf(double sand_input, double clay_input){
+ int i,index;
+ double temp, hf;
+ double silt_input=0.0; //Rawls et al (1990)
+ //do not have silt input
+ // set up mark index for inside/outside polygon check
+ double mark[POLYGON_DIMENSION]={0.0};
+ //printf("in prct2hf(), cm\n");
+ //setup the 3Dvectors and initialize them
+ struct vector cls[POLYGON_DIMENSION] = {0.0};
+ //In case silt is not == 0.0, fill up explicitly
+ for(i=0;i<POLYGON_DIMENSION;i++){
+ cls[0].sand=0.0;
+ cls[0].clay=0.0;
+ cls[0].silt=0.0;
+ }
+ cls[0].sand=0.0;
+ cls[0].clay=100.0;
+ cls[1].sand=0.0;
+ cls[1].clay=54.0;
+ cls[2].sand=17.0;
+ cls[2].clay=55.0;
+ //Get started
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ hf=175.0;
+ index=1;
+ //printf("hf=175.0\n");
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=0.0;
+ cls[0].clay=100.0;
+ cls[1].sand=17.0;
+ cls[1].clay=55.0;
+ cls[2].sand=30.0;
+ cls[2].clay=60.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=175.0\n");
+ hf=175.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=0.0;
+ cls[0].clay=100.0;
+ cls[1].sand=34.0;
+ cls[1].clay=66.0;
+ cls[2].sand=30.0;
+ cls[2].clay=60.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=175.0\n");
+ hf=175.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=100.0;
+ cls[0].clay=0.0;
+ cls[1].sand=72.0;
+ cls[1].clay=0.0;
+ cls[2].sand=65.0;
+ cls[2].clay=15.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=5.0\n");
+ hf=5.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=100.0;
+ cls[0].clay=0.0;
+ cls[1].sand=65.0;
+ cls[1].clay=30.0;
+ cls[2].sand=65.0;
+ cls[2].clay=15.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=5.0\n");
+ hf=5.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=100.0;
+ cls[0].clay=0.0;
+ cls[1].sand=65.0;
+ cls[1].clay=30.0;
+ cls[2].sand=67.0;
+ cls[2].clay=33.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=5.0\n");
+ hf=5.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=58.0;
+ cls[0].clay=42.0;
+ cls[1].sand=65.0;
+ cls[1].clay=30.0;
+ cls[2].sand=67.0;
+ cls[2].clay=33.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=15.0\n");
+ hf=15.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=58.0;
+ cls[0].clay=42.0;
+ cls[1].sand=65.0;
+ cls[1].clay=30.0;
+ cls[2].sand=30.0;
+ cls[2].clay=0.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=15.0\n");
+ hf=15.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=65.0;
+ cls[0].clay=15.0;
+ cls[1].sand=65.0;
+ cls[1].clay=30.0;
+ cls[2].sand=30.0;
+ cls[2].clay=0.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=15.0\n");
+ hf=15.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=65.0;
+ cls[0].clay=15.0;
+ cls[1].sand=72.0;
+ cls[1].clay=0.0;
+ cls[2].sand=30.0;
+ cls[2].clay=0.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=15.0\n");
+ hf=15.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=16.0;
+ cls[0].clay=0.0;
+ cls[1].sand=20.0;
+ cls[1].clay=14.0;
+ cls[2].sand=30.0;
+ cls[2].clay=0.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=25.0\n");
+ hf=25.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=28.0;
+ cls[0].clay=24.0;
+ cls[1].sand=20.0;
+ cls[1].clay=14.0;
+ cls[2].sand=30.0;
+ cls[2].clay=0.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=25.0\n");
+ hf=25.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=28.0;
+ cls[0].clay=24.0;
+ cls[1].sand=58.0;
+ cls[1].clay=42.0;
+ cls[2].sand=30.0;
+ cls[2].clay=0.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=25.0\n");
+ hf=25.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=28.0;
+ cls[0].clay=24.0;
+ cls[1].sand=58.0;
+ cls[1].clay=42.0;
+ cls[2].sand=55.0;
+ cls[2].clay=45.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=25.0\n");
+ hf=25.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=55.0;
+ cls[0].clay=45.0;
+ cls[1].sand=50.0;
+ cls[1].clay=50.0;
+ cls[2].sand=35.0;
+ cls[2].clay=42.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=35.0\n");
+ hf=35.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=55.0;
+ cls[0].clay=45.0;
+ cls[1].sand=28.0;
+ cls[1].clay=24.0;
+ cls[2].sand=35.0;
+ cls[2].clay=42.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=35.0\n");
+ hf=35.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=20.0;
+ cls[0].clay=28.0;
+ cls[1].sand=28.0;
+ cls[1].clay=24.0;
+ cls[2].sand=35.0;
+ cls[2].clay=42.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=35.0\n");
+ hf=35.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=20.0;
+ cls[0].clay=28.0;
+ cls[1].sand=28.0;
+ cls[1].clay=24.0;
+ cls[2].sand=20.0;
+ cls[2].clay=14.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=35.0\n");
+ hf=35.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=20.0;
+ cls[0].clay=28.0;
+ cls[1].sand=10.0;
+ cls[1].clay=14.0;
+ cls[2].sand=20.0;
+ cls[2].clay=14.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=35.0\n");
+ hf=35.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=16.0;
+ cls[0].clay=0.0;
+ cls[1].sand=10.0;
+ cls[1].clay=14.0;
+ cls[2].sand=20.0;
+ cls[2].clay=14.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=35.0\n");
+ hf=35.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=16.0;
+ cls[0].clay=0.0;
+ cls[1].sand=10.0;
+ cls[1].clay=14.0;
+ cls[2].sand=7.0;
+ cls[2].clay=3.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=35.0\n");
+ hf=35.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=16.0;
+ cls[0].clay=0.0;
+ cls[1].sand=0.0;
+ cls[1].clay=0.0;
+ cls[2].sand=7.0;
+ cls[2].clay=3.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=35.0\n");
+ hf=35.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=3.0;
+ cls[0].clay=3.0;
+ cls[1].sand=0.0;
+ cls[1].clay=0.0;
+ cls[2].sand=7.0;
+ cls[2].clay=3.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=35.0\n");
+ hf=35.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=3.0;
+ cls[0].clay=3.0;
+ cls[1].sand=0.0;
+ cls[1].clay=0.0;
+ cls[2].sand=0.0;
+ cls[2].clay=9.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=35.0\n");
+ hf=35.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=3.0;
+ cls[0].clay=3.0;
+ cls[1].sand=0.0;
+ cls[1].clay=28.0;
+ cls[2].sand=0.0;
+ cls[2].clay=9.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=50.0\n");
+ hf=50.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=3.0;
+ cls[0].clay=3.0;
+ cls[1].sand=0.0;
+ cls[1].clay=28.0;
+ cls[2].sand=7.0;
+ cls[2].clay=3.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=50.0\n");
+ hf=50.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=10.0;
+ cls[0].clay=14.0;
+ cls[1].sand=0.0;
+ cls[1].clay=28.0;
+ cls[2].sand=7.0;
+ cls[2].clay=3.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=50.0\n");
+ hf=50.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=10.0;
+ cls[0].clay=14.0;
+ cls[1].sand=0.0;
+ cls[1].clay=28.0;
+ cls[2].sand=10.0;
+ cls[2].clay=27.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=50.0\n");
+ hf=50.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=10.0;
+ cls[0].clay=14.0;
+ cls[1].sand=20.0;
+ cls[1].clay=28.0;
+ cls[2].sand=10.0;
+ cls[2].clay=27.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=50.0\n");
+ hf=50.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=14.0;
+ cls[0].clay=30.0;
+ cls[1].sand=20.0;
+ cls[1].clay=28.0;
+ cls[2].sand=10.0;
+ cls[2].clay=27.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=50.0\n");
+ hf=50.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=14.0;
+ cls[0].clay=30.0;
+ cls[1].sand=20.0;
+ cls[1].clay=28.0;
+ cls[2].sand=16.0;
+ cls[2].clay=35.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=50.0\n");
+ hf=50.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=23.0;
+ cls[0].clay=42.0;
+ cls[1].sand=20.0;
+ cls[1].clay=28.0;
+ cls[2].sand=16.0;
+ cls[2].clay=35.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=50.0\n");
+ hf=50.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=23.0;
+ cls[0].clay=42.0;
+ cls[1].sand=20.0;
+ cls[1].clay=28.0;
+ cls[2].sand=35.0;
+ cls[2].clay=42.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=50.0\n");
+ hf=50.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=23.0;
+ cls[0].clay=42.0;
+ cls[1].sand=36.0;
+ cls[1].clay=46.0;
+ cls[2].sand=35.0;
+ cls[2].clay=42.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=50.0\n");
+ hf=50.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=50.0;
+ cls[0].clay=50.0;
+ cls[1].sand=36.0;
+ cls[1].clay=46.0;
+ cls[2].sand=35.0;
+ cls[2].clay=42.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=50.0\n");
+ hf=50.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=50.0;
+ cls[0].clay=50.0;
+ cls[1].sand=36.0;
+ cls[1].clay=46.0;
+ cls[2].sand=45.0;
+ cls[2].clay=55.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=50.0\n");
+ hf=50.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=38.0;
+ cls[0].clay=55.0;
+ cls[1].sand=36.0;
+ cls[1].clay=46.0;
+ cls[2].sand=45.0;
+ cls[2].clay=55.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=80.0\n");
+ hf=80.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=38.0;
+ cls[0].clay=55.0;
+ cls[1].sand=42.0;
+ cls[1].clay=58.0;
+ cls[2].sand=45.0;
+ cls[2].clay=55.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=80.0\n");
+ hf=80.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=38.0;
+ cls[0].clay=55.0;
+ cls[1].sand=36.0;
+ cls[1].clay=46.0;
+ cls[2].sand=23.0;
+ cls[2].clay=42.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=80.0\n");
+ hf=80.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=38.0;
+ cls[0].clay=55.0;
+ cls[1].sand=0.0;
+ cls[1].clay=44.0;
+ cls[2].sand=23.0;
+ cls[2].clay=42.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=80.0\n");
+ hf=80.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=16.0;
+ cls[0].clay=35.0;
+ cls[1].sand=0.0;
+ cls[1].clay=44.0;
+ cls[2].sand=23.0;
+ cls[2].clay=42.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=80.0\n");
+ hf=80.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=16.0;
+ cls[0].clay=35.0;
+ cls[1].sand=0.0;
+ cls[1].clay=44.0;
+ cls[2].sand=14.0;
+ cls[2].clay=30.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=80.0\n");
+ hf=80.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=10.0;
+ cls[0].clay=27.0;
+ cls[1].sand=0.0;
+ cls[1].clay=44.0;
+ cls[2].sand=14.0;
+ cls[2].clay=30.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=80.0\n");
+ hf=80.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=10.0;
+ cls[0].clay=27.0;
+ cls[1].sand=0.0;
+ cls[1].clay=44.0;
+ cls[2].sand=0.0;
+ cls[2].clay=28.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=80.0\n");
+ hf=80.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=17.0;
+ cls[0].clay=55.0;
+ cls[1].sand=0.0;
+ cls[1].clay=44.0;
+ cls[2].sand=0.0;
+ cls[2].clay=54.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=125.0\n");
+ hf=125.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=17.0;
+ cls[0].clay=55.0;
+ cls[1].sand=0.0;
+ cls[1].clay=44.0;
+ cls[2].sand=38.0;
+ cls[2].clay=55.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=125.0\n");
+ hf=125.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=17.0;
+ cls[0].clay=55.0;
+ cls[1].sand=30.0;
+ cls[1].clay=60.0;
+ cls[2].sand=38.0;
+ cls[2].clay=55.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=125.0\n");
+ hf=125.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=42.0;
+ cls[0].clay=58.0;
+ cls[1].sand=30.0;
+ cls[1].clay=60.0;
+ cls[2].sand=38.0;
+ cls[2].clay=55.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=125.0\n");
+ hf=125.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=42.0;
+ cls[0].clay=58.0;
+ cls[1].sand=30.0;
+ cls[1].clay=60.0;
+ cls[2].sand=34.0;
+ cls[2].clay=66.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("hf=125.0\n");
+ hf=125.0;
+ }
+ }
+ return hf;
+}
Added: grass-addons/gipe/r.soiltex2prop/prct2ksat.c
===================================================================
--- grass-addons/gipe/r.soiltex2prop/prct2ksat.c (rev 0)
+++ grass-addons/gipe/r.soiltex2prop/prct2ksat.c 2008-05-22 13:41:51 UTC (rev 31483)
@@ -0,0 +1,1038 @@
+#include<stdio.h>
+
+#define POLYGON_DIMENSION 50
+
+struct vector{
+ double sand;
+ double clay;
+ double silt;
+};
+
+
+// KSAT
+
+int prct2ksat(double sand_input, double clay_input){
+ int i,index;
+ double temp,ksat;
+ double silt_input=0.0; //Rawls et al (1990)
+ //do not have silt input
+ // set up mark index for inside/outside polygon check
+ double mark[POLYGON_DIMENSION]={0.0};
+ //printf("in prct2ksat(), cm/h\n");
+ //setup the 3Dvectors and initialize them
+ struct vector cls[POLYGON_DIMENSION] = {0.0};
+ //In case silt is not == 0.0, fill up explicitly
+ for(i=0;i<POLYGON_DIMENSION;i++){
+ cls[0].sand=0.0;
+ cls[0].clay=0.0;
+ cls[0].silt=0.0;
+ }
+ cls[0].sand=0.0;
+ cls[0].clay=100.0;
+ cls[1].sand=0.0;
+ cls[1].clay=60.0;
+ cls[2].sand=7.0;
+ cls[2].clay=56.0;
+ //Get started
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ ksat=0.0025;
+ index=1;
+ //printf("Ksat=0.0025\n");
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=0.0;
+ cls[0].clay=100.0;
+ cls[1].sand=7.0;
+ cls[1].clay=56.0;
+ cls[2].sand=30.0;
+ cls[2].clay=56.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.0025\n");
+ ksat=0.0025;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=0.0;
+ cls[0].clay=100.0;
+ cls[1].sand=30.0;
+ cls[1].clay=56.0;
+ cls[2].sand=40.0;
+ cls[2].clay=63.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.0025\n");
+ ksat=0.0025;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=100.0;
+ cls[0].clay=0.0;
+ cls[1].sand=85.0;
+ cls[1].clay=0.0;
+ cls[2].sand=90.0;
+ cls[2].clay=10.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=25.0\n");
+ ksat=25.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=85.0;
+ cls[0].clay=0.0;
+ cls[1].sand=90.0;
+ cls[1].clay=10.0;
+ cls[2].sand=80.0;
+ cls[2].clay=0.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=15.0\n");
+ ksat=15.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=80.0;
+ cls[0].clay=0.0;
+ cls[1].sand=90.0;
+ cls[1].clay=10.0;
+ cls[2].sand=85.0;
+ cls[2].clay=15.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=15.0\n");
+ ksat=15.0;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=85.0;
+ cls[0].clay=15.0;
+ cls[1].sand=80.0;
+ cls[1].clay=0.0;
+ cls[2].sand=70.0;
+ cls[2].clay=0.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=7.5\n");
+ ksat=7.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=70.0;
+ cls[0].clay=0.0;
+ cls[1].sand=85.0;
+ cls[1].clay=15.0;
+ cls[2].sand=75.0;
+ cls[2].clay=25.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=7.5\n");
+ ksat=7.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=70.0;
+ cls[0].clay=0.0;
+ cls[1].sand=75.0;
+ cls[1].clay=25.0;
+ cls[2].sand=68.0;
+ cls[2].clay=23.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=2.5\n");
+ ksat=2.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=75.0;
+ cls[0].clay=25.0;
+ cls[1].sand=68.0;
+ cls[1].clay=23.0;
+ cls[2].sand=70.0;
+ cls[2].clay=30.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=2.5\n");
+ ksat=2.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=70.0;
+ cls[0].clay=0.0;
+ cls[1].sand=35.0;
+ cls[1].clay=0.0;
+ cls[2].sand=68.0;
+ cls[2].clay=23.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=2.5\n");
+ ksat=2.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=17.0;
+ cls[0].clay=0.0;
+ cls[1].sand=35.0;
+ cls[1].clay=0.0;
+ cls[2].sand=22.0;
+ cls[2].clay=5.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=1.5\n");
+ ksat=1.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=22.0;
+ cls[0].clay=5.0;
+ cls[1].sand=35.0;
+ cls[1].clay=0.0;
+ cls[2].sand=68.0;
+ cls[2].clay=23.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=1.5\n");
+ ksat=1.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=22.0;
+ cls[0].clay=5.0;
+ cls[1].sand=68.0;
+ cls[1].clay=23.0;
+ cls[2].sand=60.0;
+ cls[2].clay=25.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=1.5\n");
+ ksat=1.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=60.0;
+ cls[0].clay=25.0;
+ cls[1].sand=68.0;
+ cls[1].clay=23.0;
+ cls[2].sand=70.0;
+ cls[2].clay=30.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=1.5\n");
+ ksat=1.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=60.0;
+ cls[0].clay=25.0;
+ cls[1].sand=65.0;
+ cls[1].clay=35.0;
+ cls[2].sand=70.0;
+ cls[2].clay=30.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=1.5\n");
+ ksat=1.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=17.0;
+ cls[0].clay=0.0;
+ cls[1].sand=10.0;
+ cls[1].clay=0.0;
+ cls[2].sand=22.0;
+ cls[2].clay=5.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.8\n");
+ ksat=0.8;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=20.0;
+ cls[0].clay=12.0;
+ cls[1].sand=10.0;
+ cls[1].clay=0.0;
+ cls[2].sand=22.0;
+ cls[2].clay=5.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.8\n");
+ ksat=0.8;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=20.0;
+ cls[0].clay=12.0;
+ cls[1].sand=42.0;
+ cls[1].clay=20.0;
+ cls[2].sand=22.0;
+ cls[2].clay=5.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.8\n");
+ ksat=0.8;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=60.0;
+ cls[0].clay=25.0;
+ cls[1].sand=42.0;
+ cls[1].clay=20.0;
+ cls[2].sand=22.0;
+ cls[2].clay=5.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.8\n");
+ ksat=0.8;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=60.0;
+ cls[0].clay=25.0;
+ cls[1].sand=42.0;
+ cls[1].clay=20.0;
+ cls[2].sand=57.0;
+ cls[2].clay=30.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.8\n");
+ ksat=0.8;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=60.0;
+ cls[0].clay=25.0;
+ cls[1].sand=65.0;
+ cls[1].clay=35.0;
+ cls[2].sand=57.0;
+ cls[2].clay=30.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.8\n");
+ ksat=0.8;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=63.0;
+ cls[0].clay=38.0;
+ cls[1].sand=65.0;
+ cls[1].clay=35.0;
+ cls[2].sand=57.0;
+ cls[2].clay=30.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.8\n");
+ ksat=0.8;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=55.0;
+ cls[0].clay=35.0;
+ cls[1].sand=60.0;
+ cls[1].clay=40.0;
+ cls[2].sand=63.0;
+ cls[2].clay=38.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.5\n");
+ ksat=0.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=55.0;
+ cls[0].clay=35.0;
+ cls[1].sand=57.0;
+ cls[1].clay=30.0;
+ cls[2].sand=63.0;
+ cls[2].clay=38.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.5\n");
+ ksat=0.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=55.0;
+ cls[0].clay=35.0;
+ cls[1].sand=57.0;
+ cls[1].clay=30.0;
+ cls[2].sand=38.0;
+ cls[2].clay=23.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.5\n");
+ ksat=0.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=42.0;
+ cls[0].clay=20.0;
+ cls[1].sand=57.0;
+ cls[1].clay=30.0;
+ cls[2].sand=38.0;
+ cls[2].clay=23.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.5\n");
+ ksat=0.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=42.0;
+ cls[0].clay=20.0;
+ cls[1].sand=23.0;
+ cls[1].clay=20.0;
+ cls[2].sand=20.0;
+ cls[2].clay=12.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.5\n");
+ ksat=0.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=7.0;
+ cls[0].clay=3.0;
+ cls[1].sand=23.0;
+ cls[1].clay=20.0;
+ cls[2].sand=20.0;
+ cls[2].clay=12.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.5\n");
+ ksat=0.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=7.0;
+ cls[0].clay=3.0;
+ cls[1].sand=10.0;
+ cls[1].clay=0.0;
+ cls[2].sand=20.0;
+ cls[2].clay=12.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.5\n");
+ ksat=0.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=7.0;
+ cls[0].clay=3.0;
+ cls[1].sand=10.0;
+ cls[1].clay=0.0;
+ cls[2].sand=0.0;
+ cls[2].clay=0.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.5\n");
+ ksat=0.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=7.0;
+ cls[0].clay=3.0;
+ cls[1].sand=0.0;
+ cls[1].clay=3.0;
+ cls[2].sand=0.0;
+ cls[2].clay=0.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.5\n");
+ ksat=0.5;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=7.0;
+ cls[0].clay=3.0;
+ cls[1].sand=0.0;
+ cls[1].clay=3.0;
+ cls[2].sand=9.0;
+ cls[2].clay=18.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.3\n");
+ ksat=0.3;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=0.0;
+ cls[0].clay=16.0;
+ cls[1].sand=0.0;
+ cls[1].clay=3.0;
+ cls[2].sand=9.0;
+ cls[2].clay=18.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.3\n");
+ ksat=0.3;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=7.0;
+ cls[0].clay=3.0;
+ cls[1].sand=23.0;
+ cls[1].clay=20.0;
+ cls[2].sand=9.0;
+ cls[2].clay=18.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.3\n");
+ ksat=0.3;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=22.0;
+ cls[0].clay=29.0;
+ cls[1].sand=23.0;
+ cls[1].clay=20.0;
+ cls[2].sand=9.0;
+ cls[2].clay=18.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.3\n");
+ ksat=0.3;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=22.0;
+ cls[0].clay=29.0;
+ cls[1].sand=23.0;
+ cls[1].clay=20.0;
+ cls[2].sand=33.0;
+ cls[2].clay=29.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.3\n");
+ ksat=0.3;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=38.0;
+ cls[0].clay=23.0;
+ cls[1].sand=23.0;
+ cls[1].clay=20.0;
+ cls[2].sand=33.0;
+ cls[2].clay=29.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.3\n");
+ ksat=0.3;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=38.0;
+ cls[0].clay=23.0;
+ cls[1].sand=50.0;
+ cls[1].clay=35.0;
+ cls[2].sand=33.0;
+ cls[2].clay=29.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.3\n");
+ ksat=0.3;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=38.0;
+ cls[0].clay=23.0;
+ cls[1].sand=50.0;
+ cls[1].clay=35.0;
+ cls[2].sand=55.0;
+ cls[2].clay=35.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.3\n");
+ ksat=0.3;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=60.0;
+ cls[0].clay=40.0;
+ cls[1].sand=50.0;
+ cls[1].clay=35.0;
+ cls[2].sand=55.0;
+ cls[2].clay=35.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.3\n");
+ ksat=0.3;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=60.0;
+ cls[0].clay=40.0;
+ cls[1].sand=50.0;
+ cls[1].clay=35.0;
+ cls[2].sand=57.0;
+ cls[2].clay=44.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.3\n");
+ ksat=0.3;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=40.0;
+ cls[0].clay=38.0;
+ cls[1].sand=50.0;
+ cls[1].clay=35.0;
+ cls[2].sand=57.0;
+ cls[2].clay=44.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.15\n");
+ ksat=0.15;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=40.0;
+ cls[0].clay=38.0;
+ cls[1].sand=55.0;
+ cls[1].clay=45.0;
+ cls[2].sand=57.0;
+ cls[2].clay=44.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.15\n");
+ ksat=0.15;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=40.0;
+ cls[0].clay=38.0;
+ cls[1].sand=50.0;
+ cls[1].clay=35.0;
+ cls[2].sand=33.0;
+ cls[2].clay=29.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.15\n");
+ ksat=0.15;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=40.0;
+ cls[0].clay=38.0;
+ cls[1].sand=30.0;
+ cls[1].clay=38.0;
+ cls[2].sand=33.0;
+ cls[2].clay=29.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.15\n");
+ ksat=0.15;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=22.0;
+ cls[0].clay=29.0;
+ cls[1].sand=30.0;
+ cls[1].clay=38.0;
+ cls[2].sand=33.0;
+ cls[2].clay=29.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.15\n");
+ ksat=0.15;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=22.0;
+ cls[0].clay=29.0;
+ cls[1].sand=30.0;
+ cls[1].clay=38.0;
+ cls[2].sand=13.0;
+ cls[2].clay=29.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.15\n");
+ ksat=0.15;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=22.0;
+ cls[0].clay=29.0;
+ cls[1].sand=9.0;
+ cls[1].clay=18.0;
+ cls[2].sand=13.0;
+ cls[2].clay=29.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.15\n");
+ ksat=0.15;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=0.0;
+ cls[0].clay=16.0;
+ cls[1].sand=9.0;
+ cls[1].clay=18.0;
+ cls[2].sand=13.0;
+ cls[2].clay=29.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.15\n");
+ ksat=0.15;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=0.0;
+ cls[0].clay=16.0;
+ cls[1].sand=0.0;
+ cls[1].clay=25.0;
+ cls[2].sand=13.0;
+ cls[2].clay=29.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.15\n");
+ ksat=0.15;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=0.0;
+ cls[0].clay=33.0;
+ cls[1].sand=0.0;
+ cls[1].clay=25.0;
+ cls[2].sand=15.0;
+ cls[2].clay=38.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.075\n");
+ ksat=0.075;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=13.0;
+ cls[0].clay=29.0;
+ cls[1].sand=0.0;
+ cls[1].clay=25.0;
+ cls[2].sand=15.0;
+ cls[2].clay=38.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.075\n");
+ ksat=0.075;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=13.0;
+ cls[0].clay=29.0;
+ cls[1].sand=30.0;
+ cls[1].clay=38.0;
+ cls[2].sand=15.0;
+ cls[2].clay=38.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.075\n");
+ ksat=0.075;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=20.0;
+ cls[0].clay=41.0;
+ cls[1].sand=30.0;
+ cls[1].clay=38.0;
+ cls[2].sand=15.0;
+ cls[2].clay=38.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.075\n");
+ ksat=0.075;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=20.0;
+ cls[0].clay=41.0;
+ cls[1].sand=30.0;
+ cls[1].clay=38.0;
+ cls[2].sand=50.0;
+ cls[2].clay=50.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.075\n");
+ ksat=0.075;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=40.0;
+ cls[0].clay=38.0;
+ cls[1].sand=30.0;
+ cls[1].clay=38.0;
+ cls[2].sand=50.0;
+ cls[2].clay=50.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.075\n");
+ ksat=0.075;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=40.0;
+ cls[0].clay=38.0;
+ cls[1].sand=55.0;
+ cls[1].clay=45.0;
+ cls[2].sand=50.0;
+ cls[2].clay=50.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.075\n");
+ ksat=0.075;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=13.0;
+ cls[0].clay=50.0;
+ cls[1].sand=20.0;
+ cls[1].clay=41.0;
+ cls[2].sand=50.0;
+ cls[2].clay=50.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.025\n");
+ ksat=0.025;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=13.0;
+ cls[0].clay=50.0;
+ cls[1].sand=18.0;
+ cls[1].clay=53.0;
+ cls[2].sand=50.0;
+ cls[2].clay=50.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.025\n");
+ ksat=0.025;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=35.0;
+ cls[0].clay=55.0;
+ cls[1].sand=18.0;
+ cls[1].clay=53.0;
+ cls[2].sand=50.0;
+ cls[2].clay=50.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.025\n");
+ ksat=0.025;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=35.0;
+ cls[0].clay=55.0;
+ cls[1].sand=43.0;
+ cls[1].clay=58.0;
+ cls[2].sand=50.0;
+ cls[2].clay=50.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.025\n");
+ ksat=0.025;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=13.0;
+ cls[0].clay=50.0;
+ cls[1].sand=20.0;
+ cls[1].clay=41.0;
+ cls[2].sand=15.0;
+ cls[2].clay=38.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.025\n");
+ ksat=0.025;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=13.0;
+ cls[0].clay=50.0;
+ cls[1].sand=0.0;
+ cls[1].clay=33.0;
+ cls[2].sand=15.0;
+ cls[2].clay=38.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.025\n");
+ ksat=0.025;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=13.0;
+ cls[0].clay=50.0;
+ cls[1].sand=0.0;
+ cls[1].clay=33.0;
+ cls[2].sand=0.0;
+ cls[2].clay=51.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.025\n");
+ ksat=0.025;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=8.0;
+ cls[0].clay=56.0;
+ cls[1].sand=0.0;
+ cls[1].clay=60.0;
+ cls[2].sand=0.0;
+ cls[2].clay=51.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.0075\n");
+ ksat=0.0075;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=8.0;
+ cls[0].clay=56.0;
+ cls[1].sand=13.0;
+ cls[1].clay=50.0;
+ cls[2].sand=0.0;
+ cls[2].clay=51.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.0075\n");
+ ksat=0.0075;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=8.0;
+ cls[0].clay=56.0;
+ cls[1].sand=13.0;
+ cls[1].clay=50.0;
+ cls[2].sand=18.0;
+ cls[2].clay=53.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.0075\n");
+ ksat=0.0075;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=8.0;
+ cls[0].clay=56.0;
+ cls[1].sand=30.0;
+ cls[1].clay=56.0;
+ cls[2].sand=18.0;
+ cls[2].clay=53.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.0075\n");
+ ksat=0.0075;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=35.0;
+ cls[0].clay=55.0;
+ cls[1].sand=30.0;
+ cls[1].clay=56.0;
+ cls[2].sand=18.0;
+ cls[2].clay=53.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.0075\n");
+ ksat=0.0075;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=35.0;
+ cls[0].clay=55.0;
+ cls[1].sand=30.0;
+ cls[1].clay=56.0;
+ cls[2].sand=43.0;
+ cls[2].clay=58.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.0075\n");
+ ksat=0.0075;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=40.0;
+ cls[0].clay=63.0;
+ cls[1].sand=30.0;
+ cls[1].clay=56.0;
+ cls[2].sand=43.0;
+ cls[2].clay=58.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Ksat=0.0075\n");
+ ksat=0.0075;
+ }
+ }
+ return ksat;
+}
Added: grass-addons/gipe/r.soiltex2prop/prct2porosity.c
===================================================================
--- grass-addons/gipe/r.soiltex2prop/prct2porosity.c (rev 0)
+++ grass-addons/gipe/r.soiltex2prop/prct2porosity.c 2008-05-22 13:41:51 UTC (rev 31483)
@@ -0,0 +1,244 @@
+#include<stdio.h>
+
+#define POLYGON_DIMENSION 50
+
+struct vector{
+ double sand;
+ double clay;
+ double silt;
+};
+
+int prct2porosity(double sand_input, double clay_input){
+ int i,index;
+ double temp,porosity;
+ double silt_input=0.0; //Rawls et al (1990)
+ //do not have silt input
+ // set up mark index for inside/outside polygon check
+ double mark[POLYGON_DIMENSION]={0.0};
+ //printf("in prct2poros(), Volume Fraction\n");
+ //setup the 3Dvectors and initialize them
+ struct vector cls[POLYGON_DIMENSION] = {0.0};
+ //In case silt is not == 0.0, fill up explicitly
+ for(i=0;i<POLYGON_DIMENSION;i++){
+ cls[0].sand=0.0;
+ cls[0].clay=0.0;
+ cls[0].silt=0.0;
+ }
+ cls[0].sand=0.0;
+ cls[0].clay=100.0;
+ cls[1].sand=10.0;
+ cls[1].clay=90.0;
+ cls[2].sand=25.0;
+ cls[2].clay=75.0;
+ //Get started
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ porosity=0.575;
+ index=1;
+ //printf("Poros=0.575\n");
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=10.0;
+ cls[0].clay=0.0;
+ cls[1].sand=20.0;
+ cls[1].clay=20.0;
+ cls[2].sand=50.0;
+ cls[2].clay=0.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Poros=0.575\n");
+ porosity=0.575;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=100.0;
+ cls[0].clay=0.0;
+ cls[1].sand=50.0;
+ cls[1].clay=50.0;
+ cls[2].sand=50.0;
+ cls[2].clay=43.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Poros=0.425\n");
+ porosity=0.425;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=100.0;
+ cls[0].clay=0.0;
+ cls[1].sand=50.0;
+ cls[1].clay=43.0;
+ cls[2].sand=52.0;
+ cls[2].clay=33.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Poros=0.425\n");
+ porosity=0.425;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=100.0;
+ cls[0].clay=0.0;
+ cls[1].sand=52.0;
+ cls[1].clay=33.0;
+ cls[2].sand=57.0;
+ cls[2].clay=25.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Poros=0.425\n");
+ porosity=0.425;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=100.0;
+ cls[0].clay=0.0;
+ cls[1].sand=57.0;
+ cls[1].clay=25.0;
+ cls[2].sand=87.0;
+ cls[2].clay=0.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Poros=0.425\n");
+ porosity=0.425;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=0.0;
+ cls[0].clay=0.0;
+ cls[1].sand=25.0;
+ cls[1].clay=75.0;
+ cls[2].sand=0.0;
+ cls[2].clay=90.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Poros=0.525\n");
+ porosity=0.525;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=0.0;
+ cls[0].clay=0.0;
+ cls[1].sand=25.0;
+ cls[1].clay=75.0;
+ cls[2].sand=10.0;
+ cls[2].clay=0.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Poros=0.525\n");
+ porosity=0.525;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=10.0;
+ cls[0].clay=0.0;
+ cls[1].sand=25.0;
+ cls[1].clay=75.0;
+ cls[2].sand=20.0;
+ cls[2].clay=20.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Poros=0.525\n");
+ porosity=0.525;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=20.0;
+ cls[0].clay=20.0;
+ cls[1].sand=25.0;
+ cls[1].clay=75.0;
+ cls[2].sand=25.0;
+ cls[2].clay=55.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Poros=0.525\n");
+ porosity=0.525;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=20.0;
+ cls[0].clay=20.0;
+ cls[1].sand=25.0;
+ cls[1].clay=55.0;
+ cls[2].sand=27.0;
+ cls[2].clay=45.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Poros=0.525\n");
+ porosity=0.525;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=20.0;
+ cls[0].clay=20.0;
+ cls[1].sand=27.0;
+ cls[1].clay=45.0;
+ cls[2].sand=50.0;
+ cls[2].clay=0.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Poros=0.525\n");
+ porosity=0.525;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=50.0;
+ cls[0].clay=0.0;
+ cls[1].sand=70.0;
+ cls[1].clay=0.0;
+ cls[2].sand=37.0;
+ cls[2].clay=25.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Poros=0.525\n");
+ porosity=0.525;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=25.0;
+ cls[0].clay=75.0;
+ cls[1].sand=25.0;
+ cls[1].clay=55.0;
+ cls[2].sand=28.0;
+ cls[2].clay=61.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Poros=0.525\n");
+ porosity=0.525;
+ }
+ }
+ if (index==0){// if index not found then continue
+ cls[0].sand=25.0;
+ cls[0].clay=75.0;
+ cls[1].sand=37.0;
+ cls[1].clay=65.0;
+ cls[2].sand=28.0;
+ cls[2].clay=61.0;
+ mark[0]=point_in_triangle(sand_input,clay_input,silt_input,cls[0].sand,cls[0].clay,cls[0].silt,cls[1].sand,cls[1].clay,cls[1].silt,cls[2].sand,cls[2].clay,cls[2].silt);
+ if(mark[0]==1){
+ index=1;
+ //printf("Poros=0.525\n");
+ porosity=0.525;
+ }
+ }
+ if (index==0){// if index not found then continue
+ index=1;
+ //printf("Poros=0.475\n");
+ porosity=0.475;
+ }
+ return porosity;
+}
+
+
Added: grass-addons/gipe/r.soiltex2prop/vector_multiplication.c
===================================================================
--- grass-addons/gipe/r.soiltex2prop/vector_multiplication.c (rev 0)
+++ grass-addons/gipe/r.soiltex2prop/vector_multiplication.c 2008-05-22 13:41:51 UTC (rev 31483)
@@ -0,0 +1,85 @@
+#include<stdio.h>
+
+#define POLYGON_DIMENSION 50
+
+struct vector{
+ double sand;
+ double clay;
+ double silt;
+};
+
+double point_in_triangle(double point_x, double point_y, double point_z, double t1_x, double t1_y, double t1_z, double t2_x, double t2_y, double t2_z, double t3_x, double t3_y, double t3_z){
+ //printf("in function: sand=%5.3f clay=%5.3f silt=%5.3f\n",point_x,point_y,point_z);
+ double answer;
+ double answer1_x, answer1_y, answer1_z;
+ double answer2_x, answer2_y, answer2_z;
+ double answer3_x, answer3_y, answer3_z;
+ // Consider three points forming a trinagle from a given soil class boundary ABC
+ // Consider F an additional point in space
+ double af1, af2, af3; //Points for vector AF
+ double bf1, bf2, bf3; //Points for vector BF
+ double cf1, cf2, cf3; //Points for vector CF
+ double ab1, ab2, ab3; //Points for vector AB
+ double bc1, bc2, bc3; //Points for vector BC
+ double ca1, ca2, ca3; //Points for vector CA
+ // Create vectors AB, BC and CA
+ ab1=(t2_x-t1_x);
+ ab2=(t2_y-t1_y);
+ ab3=(t2_z-t1_z);
+ bc1=(t3_x-t2_x);
+ bc2=(t3_y-t2_y);
+ bc3=(t3_z-t2_z);
+ ca1=(t1_x-t3_x);
+ ca2=(t1_y-t3_y);
+ ca3=(t1_z-t3_z);
+ // Create vectors AF, BF and CF
+ af1=(point_x-t1_x);
+ af2=(point_y-t1_y);
+ af3=(point_z-t1_z);
+ bf1=(point_x-t2_x);
+ bf2=(point_y-t2_y);
+ bf3=(point_z-t2_z);
+ cf1=(point_x-t3_x);
+ cf2=(point_y-t3_y);
+ cf3=(point_z-t3_z);
+ // Calculate the following CrossProducts:
+ // AFxAB
+ answer1_x=(af2*ab3)-(af3*ab2);
+ answer1_y=(af3*ab1)-(af1*ab3);
+ answer1_z=(af1*ab2)-(af2*ab1);
+// printf("answer(AFxAB)= %f %f %f\n",answer1_x, answer1_y, answer1_z);
+ //BFxBC
+ answer2_x=(bf2*bc3)-(bf3*bc2);
+ answer2_y=(bf3*bc1)-(bf1*bc3);
+ answer2_z=(bf1*bc2)-(bf2*bc1);
+// printf("answer(BFxBC)= %f %f %f\n",answer2_x, answer2_y, answer2_z);
+ //CFxCA
+ answer3_x=(cf2*ca3)-(cf3*ca2);
+ answer3_y=(cf3*ca1)-(cf1*ca3);
+ answer3_z=(cf1*ca2)-(cf2*ca1);
+// printf("answer(CFxCA)= %f %f %f\n",answer3_x, answer3_y, answer3_z);
+ answer=0.0;//initialize value
+ if((int)answer1_x>=0&&(int)answer2_x>=0&&(int)answer3_x>=0){
+ answer+=1.0;
+ } else if((int)answer1_x<=0&&(int)answer2_x<=0&&(int)answer3_x<=0){
+ answer-=1.0;
+ }
+ if((int)answer1_y>=0&&(int)answer2_y>=0&&(int)answer3_y>=0){
+ answer+=1.0;
+ } else if((int)answer1_y<=0&&(int)answer2_y<=0&&(int)answer3_y<=0){
+ answer-=1.0;
+ }
+ if((int)answer1_z>=0&&(int)answer2_z>=0&&(int)answer3_z>=0){
+ answer+=1.0;
+ } else if((int)answer1_z<=0&&(int)answer2_z<=0&&(int)answer3_z<=0){
+ answer-=1.0;
+ }
+ if(answer==3||answer==-3){
+ answer=1;
+ } else {
+ answer=0;
+ }
+ return answer;
+}
+
+
More information about the grass-commit
mailing list