/*

* C Source code

* This program computes an experimental semivariogram using common

* method described by Journel and Huijbregts, "Mining Geostatistics"

* 1978.

* It takes an input text file called "dbQuery.txt" of 5 columns:

* ID EASTING NORTHING ELEVATION CONCENTRATION

* Output is a text file "svar.o" of 4 columns:

* LAG SEMIVARIANCE ST.DEVIATION

* Written by Bart Faulkner, 1998

* RS Kerr Environmental Research Center

* PO Box 1198

* Ada, OK 74820

*

*/

#include "svar.h"

int nsvardata;

int main(void)

{

int vchk=1,hchk=1;

hchk=hsvar(90.,0.9);

return 0;

}

hsvar(double cone_ang, double csize)

{

int i,j,k,l,ll,clipto,span,entno,lastrow,nelem[maxvect];

float east[maxvect],north[maxvect],elev[maxvect],bchem[maxvect],d[1000];

double lag[maxvect],gam[maxvect],std[maxvect];

double maxeast,mineast,maxnorth,minnorth,maxelev,minelev,rad_ang,arg_ang;

double diff,sumdist[maxvect],avdist[maxvect],zdiffsq[maxvect],avzdiffsq[maxvect];

double sdevgamma[maxvect],latdiff;

FILE *ifp, *ofp;

#define hcone_side <

ifp=fopen("dbQuery.txt","r");

ofp=fopen("svar.o","wb");

rad_ang=cone_ang*3.141592654/180.;

for (i=0; i<maxvect; i++) {

east[i]=0.;

north[i]=0.;

elev[i]=0.;

bchem[i]=0.;

}

for (i=0; fscanf(ifp, "%d",&entno) != EOF; ++i) {

fscanf(ifp,"%f %f %f %f\n",&east[i],&north[i],&elev[i],&bchem[i]);

lastrow=i;

}

maxeast=0.;

mineast=1.e10;

maxnorth=0.;

minnorth=1.e10;

maxelev=0.;

minelev=1.e10;

for (i=0; i<=lastrow; i++) {

for (j=0; j<=lastrow; j++) {

if (east[i] > east[j] && east[i] > maxeast) maxeast=east[i];

if (east[i] < east[j] && east[i] < mineast) mineast=east[i];

if (north[i] > north[j] && north[i] > maxnorth) maxnorth=north[i];

if (north[i] < north[j] && north[i] < minnorth) minnorth=north[i];

if (elev[i] > elev[j] && elev[i] > maxelev) maxelev=elev[i];

if (elev[i] < elev[j] && elev[i] < minelev) minelev=elev[i];

}

}

span= (int) (maxelev-minelev)/csize;

for (i=0; i<span; i++) d[i]= (double) csize*i;

/* Calculate semivariogram: */

printf("calculating semivariogram...");

for (i=0; i<lastrow; i++) {

for (j=1; j<=lastrow; j++) {

latdiff=sqrt(pow((east[i]-east[j-1]),2.)+pow((north[i]-north[j-1]),2.));

arg_ang=atan(latdiff/fabs(elev[i]-elev[j-1]));

if (fabs(arg_ang) hcone_side fabs(rad_ang)) {

diff=sqrt(pow((east[i]-east[j-1]),2.)+pow((north[i]-north[j-1]),2.)+

pow(elev[i]-elev[j-1],2.));

for (k=1; k<span; k++) {

if (diff > d[k-1] && diff < d[k]) {

nelem[k]+=1;

sumdist[k]=sumdist[k]+diff;

avdist[k]=sumdist[k]/nelem[k];

zdiffsq[k]=zdiffsq[k]+0.5*pow((bchem[i]-bchem[j]),2.);

avzdiffsq[k]=zdiffsq[k]/nelem[k];

sdevgamma[k]=2.*sqrt((2.*pow(avzdiffsq[k],2.)/nelem[k]));

}

}

}

}

}

ll=-1;

for (l=0; l<span; l++) {

if (avdist[l] != 0. && avzdiffsq[l] != 0.) {

ll=ll+1;

lag[ll]=avdist[l];

gam[ll]=avzdiffsq[l];

std[ll]=sdevgamma[l];

clipto=ll;

}

}

for (ll=0; ll<clipto; ll++) {

fprintf(ofp,"%f %f %f\n",lag[ll],gam[ll],std[ll]);

}

nsvardata=clipto;

printf("\n%d\n",nsvardata);

return 0;

}