From 07bcff1355fc7cdf938c782618219cf33c1295b3 Mon Sep 17 00:00:00 2001 From: l2r <80277078+l2r-git@users.noreply.github.com> Date: Sat, 21 Sep 2024 10:01:54 -0700 Subject: [PATCH 1/5] Simplify parameters. Simplify initial use. Extend to any number of layers. --- src/su/main/statics/sunszonecsv.c | 517 +++++++++++++++++++----------- 1 file changed, 335 insertions(+), 182 deletions(-) diff --git a/src/su/main/statics/sunszonecsv.c b/src/su/main/statics/sunszonecsv.c index 2bce1b0e..088e4f7c 100644 --- a/src/su/main/statics/sunszonecsv.c +++ b/src/su/main/statics/sunszonecsv.c @@ -1,7 +1,7 @@ -/* Copyright (c) Colorado School of Mines, 2022.*/ +/* Copyright (c) Colorado School of Mines, 2024.*/ /* All rights reserved. */ -/* SUNSZONECSV: $Revision: 1.01 $ ; $Date: 2024/06/25 00:00:01 $ */ +/* SUNSZONECSV: $Revision: 1.02 $ ; $Date: 2024/09/06 00:00:01 $ */ #include "su.h" #include "segy.h" @@ -19,32 +19,25 @@ char *sdoc[] = { " ", " Parameters: ", " ", -" qin= Input file containing q-records. This parameter is required. ", +" qin= Input file containing q-records. If this parameter is ", +" not specified, all parameters down to thickadd are ignored. ", +" And then both thickadd and veladd must be specified. ", " ", -" qx=sx Name of X coordinate at points in the q-file. ", +" xys=sx,sy Names of 2 coordinate values in the q-file. ", +" Note: Default names match an example simple statics spreadsheet. ", " ", -" qy=sy Name of Y coordinate at points in the q-file. ", +" thicks=zw_true,zb_true Names of thickness values in the q-file. ", +" You can specify as many names as exist in the q-file but you ", +" do not have to specify all of them. ", +" Note: Default names match the example simple statics spreadsheet. ", " ", -" qvw=vw Name of top layer velocity at points in the q-file. ", -" ", -" qzw=zw_true Name of top layer thickness at points in the q-file. ", -" ", -" qvb=vb Name of second layer velocity at points in the q-file. ", -" ", -" qzb=zb_true Name of second layer thickness at points in the q-file. ", -" =unknown Second layer has unknown thickness. ", -" *** Note: Once unknown is specified (or defaults) for a thickness, ", -" parameters for deeper layers are not used. ", -" ", -" qvc=vc Name of third layer velocity at points in the q-file. ", -" ", -" qzc=zc_true Name of third layer thickness at points in the q-file. ", -" =unknown Third layer has unknown thickness. ", -" ", -" qvd=vd Name of fourth layer velocity at points in the q-file. ", +" vels=vw,vb,vc Names of velocity values in the q-file. ", +" This list must have one more name than the thicks list. ", +" Note: Default names match the example simple statics spreadsheet. ", " ", " locs=0 Use near surface model from shot and receiver ends of traces.", " =1 Use near surface model from shot end only. ", +" Note: If qin has widely spaced points, usually best to use 1 here. ", " =-1 Use near surface model from receiver end only. ", " Note: The model is assumed to be flat at both ends, and raypaths ", " are assumed to be straight with a linear change in layer ", @@ -58,15 +51,52 @@ char *sdoc[] = { " for both ends of the raypaths. Option -1 simply uses the ", " velocities and thicknesses of the receiver location for ", " both ends of the raypaths. ", +" *** Due to the approximations used herein, traveltimes for layers", +" ABOVE the shot hole depth will be too large. Traveltimes for ", +" layers CONTAINING or BELOW the shot hole depth should have ", +" similar accuracy as for surface shots. ", +" ", +" thickadd=0.0 List of values to add to thicknesses at both ends. ", +" If qin file is specified, default is 0.0 for all layers. ", +" If the result for a layer is less than 0.0, it is reset ", +" to 0.0 at the shot or receiver end (or both ends). ", +" Values listed for more than qin layers will be ignored. ", +" Note: If qin file is not specified, this list must be specified ", +" and are the thicknesses at both shot and receiver ends. ", +" ", +" veladd=0.0 List of values to add to velocities at both ends. ", +" If qin file is specified, default is 0.0 for all layers. ", +" If the result for a layer is less than 1.0, it is reset ", +" to 1.0 at the shot or receiver end (or both ends). ", +" Values listed for more than qin layers+1 will be ignored. ", +" Note: If qin file is not specified, this list must be specified ", +" and are the velocities at both shot and receiver ends. ", +" Must be one more value listed than in the thickadd list. ", +" ", +" offsub=0.0 Value to subtract from all distances from shot to receiver. ", +" Note: This program computes the offset distance using sx,sy,gx,gy ", +" (it does not use offset key value at all). The value here is ", +" subtracted from that distance. For 2d surveys this can be ", +" used to adjust the distance from the centre of the geophone ", +" array to the geophone nearest the shot. ", +" ", +" laymax= Maximum layers to use. When qin file is input, the default ", +" is number of values in vels list. When qin is not input, ", +" default is number of values in veladd list. ", +" =1 Set to top layer (direct arrival layer). ", +" =2 Set to second layer (first refracted arrival layer). ", +" =n Set to nth layer (n-1th refracted arrival layer). ", +" Note: This parameter exists simply to allow you to fully specify ", +" the thicks, vels, thickadd, veladd lists without having to ", +" change them to test results using fewer layer(s). ", " ", " tpath=0 Travel path shift value. ", " The default (0) is to set the shift to the minimum time of ", " the raypaths computed for all layers and the direct arrival ", " path (also known as the first arrival or first break time). ", -" =1 Set shift to top-of-first layer direct arrival time. ", -" =2 Set shift to top-of-second layer refraction arrival time. ", -" =3 Set shift to top-of-third layer refraction arrival time. ", -" =4 Set shift to top-of-fourth layer refraction arrival time. ", +" =1 Set to top layer (direct arrival time). ", +" =2 Set to second layer (first refracted arrival time). ", +" =n Set to nth layer (n-1th refracted arrival time). ", " ", " tapp=1 Apply the total time shift determined by the tpath option ", " and also set the sstat,gstat,tstat,sec key values. ", @@ -180,8 +210,10 @@ int main(int argc, char **argv) { cwp_String *pname = NULL; cwp_String *ndims = NULL; + cwp_String *temps = NULL; int numpname = 9; int nlayer = 1; + int laymax = 3; double *pindepa = NULL; int numdind = 0; @@ -231,20 +263,29 @@ int main(int argc, char **argv) { float *oneshift=NULL; float sampi=0.0; double tshift = 0.0; - double ptm[4]; /* Note current max nlayer is 3. But +1 needed here. */ - double stm[4]; - double gtm[4]; - double eps[4]; - double zps[4]; - double vls[4]; - double epg[4]; - double zpg[4]; - double vlg[4]; + double *ptm=NULL; + double *stm=NULL; + double *gtm=NULL; + double *eps=NULL; + double *zps=NULL; + double *vls=NULL; + double *epg=NULL; + double *zpg=NULL; + double *vlg=NULL; + double *tad=NULL; + double *vad=NULL; + + int nadthicks=0; + int nadvels=0; + double *thickadd=NULL; + double *veladd=NULL; double shole = 0.0; + double ehole = 0.0; + double offsub=0.0; cwp_String *iname = NULL; int numiname = 9; - int locn[9]; + int *locn=NULL; /* hook up getpar */ initargs(argc, argv); @@ -265,9 +306,6 @@ int main(int argc, char **argv) { tree_nt[i] = 1; } - if(!getparint("locs",&locs)) locs = 0; - if(locs<-1 || locs>1) err("**** Error: locs parameter out-of-range."); - if(!getparint("apply",&iapply)) iapply = -1; if(iapply!=-1 && iapply!=1) err("**** Error: apply parameter out-of-range."); @@ -292,101 +330,177 @@ int main(int argc, char **argv) { /*-------------------------------------------------------------------------- */ - getparstring("qin", &Pname); + locs = -2; /* set this flag to -2 when no qin file is input. */ + + if(countparval("qin") > 0) { - fpP = fopen(Pname, "r"); - if(fpP==NULL) err("error: input Q-file did not open correctly."); + getparstring("qin", &Pname); + + fpP = fopen(Pname, "r"); + if(fpP==NULL) err("error: input Q-file did not open correctly."); /* Set input numpname,pname to store everything in q-file. */ /* numpname>0 is a flag to ONLY store values if they are on pname list. */ /* numpname<1 is a flag to NOT store values if they are on pname list. */ - pname = ealloc1(9,sizeof(cwp_String *)); + if((countparval("thicks")>0 && countparval("vels") < 1) || + (countparval("thicks")<1 && countparval("vels") > 0)) + err("error: vels list and thicks list must both default or both be specified."); - if(countparval("qx") > 0) { - getparstring("qx",pname); - } - else { - pname[0] = ealloc1(2,1); - strcpy(pname[0],"sx"); - } - - if(countparval("qy") > 0) { - getparstring("qy",pname+1); - } - else { - pname[1] = ealloc1(2,1); - strcpy(pname[1],"sy"); - } - - if(countparval("qvw") > 0) { - getparstring("qvw",pname+2); - } - else { - pname[2] = ealloc1(2,1); - strcpy(pname[2],"vw"); - } - - if(countparval("qzw") > 0) { - getparstring("qzw",pname+3); - } - else { - pname[3] = ealloc1(7,1); - strcpy(pname[3],"zw_true"); - } - - if(countparval("qvb") > 0) { - getparstring("qvb",pname+4); - } - else { - pname[4] = ealloc1(2,1); - strcpy(pname[4],"vb"); - } - - if(countparval("qzb") > 0) { - getparstring("qzb",pname+5); - } - else { - pname[5] = ealloc1(7,1); - strcpy(pname[5],"zb_true"); - } - - if(countparval("qvc") > 0) { - getparstring("qvc",pname+6); - } + if(countparval("thicks")==0) { + nlayer = 2; + } + else { + nlayer = countparval("thicks"); + if(countparval("vels") != nlayer+1) err("error: vels list must have one more than thicks list."); + } + + if(!getparint("laymax",&laymax)) laymax = nlayer+1; + if(laymax<1 || laymax>nlayer+1) err("**** Error: laymax parameter out-of-range."); + nlayer = laymax-1; + + pname = ealloc1(nlayer+3+nlayer,sizeof(cwp_String *)); /* 2 for XYs, and 1 for extra vel */ + + temps = ealloc1(nlayer+3,sizeof(cwp_String *)); + + if(countparval("xys") > 0) { + if(countparval("xys") != 2) err("error: xys list must have 2 names."); + getparstringarray("xys", temps); + for(i=0; i<2; i++) { + pname[i] = ealloc1(strlen(temps[i]),1); + strcpy(pname[i],temps[i]); + } + } + else { + pname[0] = ealloc1(2,1); + strcpy(pname[0],"sx"); + pname[1] = ealloc1(2,1); + strcpy(pname[1],"sy"); + } + + if(countparval("vels") > 0) { + getparstringarray("vels", temps); + for(i=2; i 0) { + getparstringarray("thicks", temps); + for(i=nlayer+3;i1) err("**** Error: locs parameter out-of-range."); + + if(countparval("thickadd")>0) { + nadthicks = countparval("thickadd"); + thickadd = ealloc1double(nadthicks); + getpardouble("thickadd",thickadd); + if(nadthicks>nlayer) nadthicks = nlayer; + } + + if(countparval("veladd")>0) { + nadvels = countparval("veladd"); + veladd = ealloc1double(nadvels); + getpardouble("veladd",veladd); + if(nadvels>nlayer+1) nadvels = nlayer+1; + } + + } /* end of if(countparval("qin") > 0) { */ else { - pname[6] = ealloc1(2,1); - strcpy(pname[6],"vc"); - } - - if(countparval("qzc") > 0) { - getparstring("qzc",pname+7); + + nadthicks = countparval("thickadd"); + thickadd = ealloc1double(nadthicks); + getpardouble("thickadd",thickadd); + if(nadthicks<1) err("**** Error: When no qin file, the thickadd list must contain some values."); + + nlayer = nadthicks; + if(!getparint("laymax",&laymax)) laymax = nlayer+1; + if(laymax<1 || laymax>nlayer+1) err("**** Error: laymax parameter out-of-range."); + nlayer = laymax-1; + + if(nadthicks>nlayer) nadthicks = nlayer; + for(i=0; inlayer+1) nadvels = nlayer+1; + if(nadvels<=nadthicks) err("**** Error: when no qin file input, veladd list must have one more than thickadd."); + for(i=0; i 0) { - getparstring("qvd",pname+8); + + if(locs==-2) { + for(i=0; inlayer+1) err("**** Error: tpath parameter out-of-range."); @@ -396,62 +510,69 @@ int main(int argc, char **argv) { if(!getpardouble("tbulk",&tbulk)) tbulk = 200.0; -/* Note: positive numpname does not mean values come back in pname order. */ -/* Note: and returned numpname is something else, so set numiname. */ - - iname = ealloc1(numpname,sizeof(cwp_String *)); - numiname= numpname; - for(i=0; i100) - err("getviaqfile error: record %d (wrong comma count, damaged, non-numbers, ...)",errwarn-99); - else if(errwarn>0) err("getviaqfile error: unrecognized error code %d",errwarn); +/* Read-in qin file? */ + if(locs!=-2) { - checkpars(); +/* Note: positive numpname does not mean values come back in pname order. */ +/* Note: and returned numpname is something else, so set numiname. */ - if(ifixd==0) err("error: input with varying number of tuples is not allowed."); + iname = ealloc1(numpname,sizeof(cwp_String *)); + numiname= numpname; + for(i=0; i100) + err("getviaqfile error: record %d (wrong comma count, damaged, non-numbers, ...)",errwarn-99); + else if(errwarn>0) err("getviaqfile error: unrecognized error code %d",errwarn); + + if(ifixd==0) err("error: input with varying number of tuples is not allowed."); + + locn = ealloc1int(numiname); + + for(i=0; i=0) { + if(locs==0 || locs==1) { if(sdadd==1) sdist2 = sdist + sqrt(near_dist_s); cycle_for_near(tree_nodes,tree_dl,tree_nt,tree_numd, extent_min, extent_max, tars, @@ -511,14 +638,20 @@ int main(int argc, char **argv) { /* Note that nlayer is the number of thicknesses. There is an extra velocity. */ /* The shot point elevation eps[0] is from trace header. */ - for(i=0; i0.0) { + ehole = eps[0] - shole; for(i=0; i=zps[i]) { shole -= zps[i]; @@ -569,7 +700,29 @@ int main(int argc, char **argv) { break; } } - } + +/* Computing times between path ends for layers ABOVE the hole bottom is */ +/* complicated. Here, just set the layer top elevation to the hole bottom if */ +/* higher. So layers above that should produce times that are too slow. */ +/* So the first arrival time (tpath=0) should be usually be correct, and the */ +/* time for the layer CONTAINING the hole bottom should be correct. And also, */ +/* if the hole bottom is within the top layer, the direct arrival time will */ +/* be computed using true distance from the hole bottom to receiver elevation.*/ + + for(i=0; iehole) eps[i] = ehole; + } + } /* end of if(shole>0.0) { */ + +/* Compute the time between the paths ends. This is distance between shot and */ +/* receiver considering differences in X,Y, and top-of-layer depths divided */ +/* by the average velocity of each top-of-layer at the shot and the receiver. */ +/* (Averaging these two is same as linear velocity change from one to other). */ + + poff = sqrt((tars[0]-targ[0])*(tars[0]-targ[0]) + (tars[1]-targ[1])*(tars[1]-targ[1])) - offsub; + if(poff<0.0) poff = 0.0; + + for(i=0; i Date: Sat, 21 Sep 2024 10:05:16 -0700 Subject: [PATCH 2/5] Rework tests for new parameterization of sunszonecsv. --- src/demos/Geom3D/Sunszonecsv/README | 10 +++++--- .../Geom3D/Sunszonecsv/Sunszonecsv_test1 | 9 ++++--- .../Geom3D/Sunszonecsv/Sunszonecsv_test2 | 6 ++--- .../Geom3D/Sunszonecsv/Sunszonecsv_test3 | 6 ++--- .../Geom3D/Sunszonecsv/Sunszonecsv_test4 | 11 ++++---- .../Geom3D/Sunszonecsv/Sunszonecsv_test5 | 23 +++++++++++++++++ .../Geom3D/Sunszonecsv/Sunszonecsv_test6 | 25 +++++++++++++++++++ src/demos/Geom3D/Sunszonecsv/l2scrookred3.csv | 6 ++--- src/demos/Geom3D/Sunszonecsv/l2scrookred4.csv | 8 +++--- 9 files changed, 79 insertions(+), 25 deletions(-) create mode 100644 src/demos/Geom3D/Sunszonecsv/Sunszonecsv_test5 create mode 100644 src/demos/Geom3D/Sunszonecsv/Sunszonecsv_test6 diff --git a/src/demos/Geom3D/Sunszonecsv/README b/src/demos/Geom3D/Sunszonecsv/README index 7eaf8f37..46341161 100644 --- a/src/demos/Geom3D/Sunszonecsv/README +++ b/src/demos/Geom3D/Sunszonecsv/README @@ -4,13 +4,13 @@ SUNSZONECSV - Near Surface Zone Extract and Shift. The scripts in this directory are designed to test the parameters. They are not realistic tests. Notice that the input fake seismic is just 1 flat horizon. - For 4 refraction layers, actual seismic should contain up to 4 refraction - wavelets which cross-over each other at various offset distances, times. + For 3 refraction layers, actual seismic should contain up to 3 wavelets + which cross-over each other at various offset distances, times. The SUNSZONECSV program can be used either to help check simple near surface models, or for checking geometry update (are XYs correct, have traces been assigned to the correct points/stations). - Sometime near the end of 2024 I will produce a YouTube video about all this. + Sometime near the end of 2024 I will produce a YouTube video about this. ------------------------------------------------------------------------------ I advise running the scripts in the following order: @@ -26,6 +26,10 @@ sh ./Sunszonecsv_test3 sh ./Sunszonecsv_test4 +sh ./Sunszonecsv_test5 + +sh ./Sunszonecsv_test6 + --- --- diff --git a/src/demos/Geom3D/Sunszonecsv/Sunszonecsv_test1 b/src/demos/Geom3D/Sunszonecsv/Sunszonecsv_test1 index a042ee11..23c74f5c 100644 --- a/src/demos/Geom3D/Sunszonecsv/Sunszonecsv_test1 +++ b/src/demos/Geom3D/Sunszonecsv/Sunszonecsv_test1 @@ -1,10 +1,11 @@ #!/bin/sh -# Sushiftcsv_test1 - for program SUSHIFTCSV -# Author: Andre Latour, Aug 2023 +# Sunszonecsv_test1 - for program SUNSZONECSV +# Author: Andre Latour, Sept 2024 # echo "----------------------------------------------------------------------------" echo "--- Input fake traces, apply shifts. " echo "--- Basic test of computing/applying shifts (reverse and then forward). " + echo "--- Note the forward sunszonecsv also tests defaults for some parameters. " echo "----------------------------------------------------------------------------" # suxwigb fake11.su locs=0 tpath=0 tapp=1 apply=1 \ qin=l2scrookred.csv # - sugethw fake11.txt # suxwigb fake12.su locs=0 tpath=0 tapp=1 apply=-1 \ + sunszonecsv fake12.su \ qin=l2scrookred.csv # suxwigb fake22.su locs=0 tpath=0 tapp=1 apply=1 \ qin=l2scrookred2.csv # - sugethw fake22.txt # suxwigb fake32.su locs=0 tpath=0 tapp=1 apply=1 \ qin=l2scrookred3.csv # - sugethw fake32.txt # suxwigb fake41.su key1=gelev,sdepth key2=selev,sdepth key3=selev,sdepth a=0,120 # sunszonecsv fake42.su locs=0 tpath=0 tapp=1 apply=1 \ - qx=skkx qy=sandrey qvw=vaw qzw=zbbw_true qvb=vnnnb qzb=zbbw_true \ - qvc=vwwwc qzc=zccc_true qvd=vddd \ + xys=skkx,sandrey thicks=zbbw_true,zkkb_true,yyzc_true \ + vels=vaw,vnqb,vwwwc,vddd \ qin=l2scrookred4.csv # - sugethw fake42.txt # suxwigb fake51.su tpath=0 tapp=1 apply=1 \ + thickadd=15,100,200 veladd=610,1200,2200,3000 offsub=100 +# + sugethw fake51.txt +# + suxwigb fake52.su tpath=0 tapp=1 apply=-1 \ + thickadd=15,100,200 veladd=610,1200,2200,3000 offsub=100 +# + suxwigb fake61.su tpath=0 tapp=1 apply=1 \ + qin=l2scrookred.csv \ + thickadd=0,100,0 veladd=0,0,1000,0 offsub=100 +# + sugethw fake61.txt +# + suxwigb fake62.su tpath=0 tapp=1 apply=-1 \ + qin=l2scrookred.csv \ + thickadd=0,100,0 veladd=0,0,1000,0 offsub=100 +# + suxwigb Date: Sat, 9 Nov 2024 11:23:41 -0800 Subject: [PATCH 3/5] Create README --- src/demos/Geom3D/Sunearqcsv/README | 1 + 1 file changed, 1 insertion(+) create mode 100644 src/demos/Geom3D/Sunearqcsv/README diff --git a/src/demos/Geom3D/Sunearqcsv/README b/src/demos/Geom3D/Sunearqcsv/README new file mode 100644 index 00000000..3b62dc18 --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/README @@ -0,0 +1 @@ +Start From bbe653c746dd2776882d246d4005e714a1dc66af Mon Sep 17 00:00:00 2001 From: l2r <80277078+l2r-git@users.noreply.github.com> Date: Sat, 9 Nov 2024 11:26:02 -0800 Subject: [PATCH 4/5] Tests for new sunearqcsv program. This program copies nearest values from one q-file into another. --- src/demos/Geom3D/Sunearqcsv/Clean.sh | 6 ++++ src/demos/Geom3D/Sunearqcsv/README | 31 +++++++++++++++++- src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test1 | 30 ++++++++++++++++++ src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test2 | 25 +++++++++++++++ src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test3 | 33 ++++++++++++++++++++ src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test4 | 23 ++++++++++++++ src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test6 | 27 ++++++++++++++++ src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test7 | 27 ++++++++++++++++ src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test8 | 32 +++++++++++++++++++ src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test9 | 29 +++++++++++++++++ src/demos/Geom3D/Sunearqcsv/base1.csv | 11 +++++++ src/demos/Geom3D/Sunearqcsv/base2.csv | 11 +++++++ src/demos/Geom3D/Sunearqcsv/fixed7.csv | 18 +++++++++++ src/demos/Geom3D/Sunearqcsv/fixedtin7.csv | 16 ++++++++++ src/demos/Geom3D/Sunearqcsv/notuple8.csv | 15 +++++++++ src/demos/Geom3D/Sunearqcsv/qx1.csv | 17 ++++++++++ src/demos/Geom3D/Sunearqcsv/qx5.csv | 21 +++++++++++++ src/demos/Geom3D/Sunearqcsv/target1.csv | 11 +++++++ src/demos/Geom3D/Sunearqcsv/target2.csv | 9 ++++++ src/demos/Geom3D/Sunearqcsv/target3.csv | 9 ++++++ src/demos/Geom3D/Sunearqcsv/vary6.csv | 17 ++++++++++ src/demos/Geom3D/Sunearqcsv/varytin6.csv | 15 +++++++++ 22 files changed, 432 insertions(+), 1 deletion(-) create mode 100644 src/demos/Geom3D/Sunearqcsv/Clean.sh create mode 100644 src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test1 create mode 100644 src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test2 create mode 100644 src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test3 create mode 100644 src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test4 create mode 100644 src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test6 create mode 100644 src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test7 create mode 100644 src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test8 create mode 100644 src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test9 create mode 100644 src/demos/Geom3D/Sunearqcsv/base1.csv create mode 100644 src/demos/Geom3D/Sunearqcsv/base2.csv create mode 100644 src/demos/Geom3D/Sunearqcsv/fixed7.csv create mode 100644 src/demos/Geom3D/Sunearqcsv/fixedtin7.csv create mode 100644 src/demos/Geom3D/Sunearqcsv/notuple8.csv create mode 100644 src/demos/Geom3D/Sunearqcsv/qx1.csv create mode 100644 src/demos/Geom3D/Sunearqcsv/qx5.csv create mode 100644 src/demos/Geom3D/Sunearqcsv/target1.csv create mode 100644 src/demos/Geom3D/Sunearqcsv/target2.csv create mode 100644 src/demos/Geom3D/Sunearqcsv/target3.csv create mode 100644 src/demos/Geom3D/Sunearqcsv/vary6.csv create mode 100644 src/demos/Geom3D/Sunearqcsv/varytin6.csv diff --git a/src/demos/Geom3D/Sunearqcsv/Clean.sh b/src/demos/Geom3D/Sunearqcsv/Clean.sh new file mode 100644 index 00000000..a6672fd7 --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/Clean.sh @@ -0,0 +1,6 @@ +#! /bin/sh +# clean up files made by tests + +rm -f qout* + +exit diff --git a/src/demos/Geom3D/Sunearqcsv/README b/src/demos/Geom3D/Sunearqcsv/README index 3b62dc18..65777a63 100644 --- a/src/demos/Geom3D/Sunearqcsv/README +++ b/src/demos/Geom3D/Sunearqcsv/README @@ -1 +1,30 @@ -Start + +SUNEARQCSV - Assign Values To Target Q-file From Nearest In Another Q-file. + +This program is similar to sunearcsv.c which +assigns values to trace header keys from nearest records in a q-file. + + +------------------------------------------------------------------------------ + I advise running the scripts in the following order. +------------------------------------------------------------------------------ + +sh ./Sunearqcsv_test1 + +sh ./Sunearqcsv_test2 + +sh ./Sunearqcsv_test3 + +sh ./Sunearqcsv_test4 + + +--- More realistic tests follow.. + +sh ./Sunearqcsv_test6 + +sh ./Sunearqcsv_test7 + +sh ./Sunearqcsv_test8 + +sh ./Sunearqcsv_test9 + diff --git a/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test1 b/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test1 new file mode 100644 index 00000000..9bad6122 --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test1 @@ -0,0 +1,30 @@ +#!/bin/sh +# Sunearqcsv_test1 - for program SUNEARQCSV +# Author: Andre Latour, Nov 2024 +# + echo "---------------------------------------------------------------------------" + echo "--- If you look at base1 you can easily see that the numbers are unique " + echo "--- and you can tell what columns and rows they come from. " + echo "--- Base2 has exactly the same numbers but the rows have been rearranged " + echo "--- and so have the columns INCLUDING THEIR NAMES. " + echo "--- The sunearqcsv program code is doing a lot of indexing to match names " + echo "--- between qin and tin q-files. Using the same options, base1 and base2 " + echo "--- should produce the same output - thus confirming that the indexing " + echo "--- has matched the same names despite the rearranged base2 rows/columns. " + echo "--- Note that target1 is the same as base1 except all numbers have an " + echo "--- extra 1 added to them (except column g which has an extra 5 added). " + echo "--- This makes it easy to see that oreps=d,c,m causes columns d,c,m " + echo "--- in target1 to be replaced with nearest from base1 or base2. " + echo "---------------------------------------------------------------------------" +# + sunearqcsv qin=base1.csv dimx=b dimy=e \ + tin=target1.csv timx=b timy=e \ + qout=qout001.csv oreps=d,c,m +# + sunearqcsv qin=base2.csv dimx=b dimy=e \ + tin=target1.csv timx=b timy=e \ + qout=qout001a.csv oreps=d,c,m +# + echo "--- Note, there should NOT be a difference next... " + diff qout001.csv qout001a.csv +# diff --git a/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test2 b/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test2 new file mode 100644 index 00000000..d5f8f0fe --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test2 @@ -0,0 +1,25 @@ +#!/bin/sh +# Sunearqcsv_test2 - for program SUNEARQCSV +# Author: Andre Latour, Nov 2024 +# + echo "---------------------------------------------------------------------------" + echo "--- Using base1 and target1 like test 1. " + echo "--- But finding nearest using just 1 dimension (name g). " + echo "--- Note that target g values end in 5. Which means they are exactly " + echo "--- half-way between the g values in base1. " + echo "--- For half-way, by default, the nearest is from higher numbered row. " + echo "--- But ordp=j option causes rows to be sorted internally by j resulting " + echo "--- in reversing this because j is high-to-low. Note columns c,d,m change. " + echo "---------------------------------------------------------------------------" +# + sunearqcsv qin=base1.csv dimx=g \ + tin=target1.csv timx=g \ + qout=qout002.csv oreps=m,d,c +# + sunearqcsv qin=base1.csv dimx=g \ + tin=target1.csv timx=g \ + qout=qout002a.csv oreps=m,d,c ordp=j +# + echo "--- Note, should be differences next (but easier to compare this way) " + diff qout002.csv qout002a.csv +# diff --git a/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test3 b/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test3 new file mode 100644 index 00000000..52764cdb --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test3 @@ -0,0 +1,33 @@ +#!/bin/sh +# Sunearqcsv_test3 - for program SUNEARQCSV +# Author: Andre Latour, Nov 2024 +# + echo "---------------------------------------------------------------------------" + echo "--- Using base1,base2 and target2. " + echo "--- Target2 is rearranged/missing rows and renamed/missing some columns. " + echo "--- Remember, when oreps name exists in target it is REPLACED in output. " + echo "--- But if oreps name is missing from target, it is ADDED to output. " + echo "---------------------------------------------------------------------------" +# + sunearqcsv qin=base1.csv dimx=f \ + tin=target2.csv timx=f \ + qout=qout003.csv oreps=e,m,d,el,c,eye +# + sunearqcsv qin=base2.csv dimx=f \ + tin=target2.csv timx=f \ + qout=qout003a.csv oreps=e,m,d,el,c,eye +# + echo " Note, there should NOT be a difference next.... " + diff qout003.csv qout003a.csv +# + sunearqcsv qin=base1.csv dimx=f \ + tin=target2.csv timx=f \ + qout=qout003b.csv oreps=c +# + sunearqcsv qin=base2.csv dimx=f \ + tin=target2.csv timx=f \ + qout=qout003c.csv oreps=c +# + echo " Note, there should NOT be a difference next.... " + diff qout003b.csv qout003c.csv +# diff --git a/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test4 b/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test4 new file mode 100644 index 00000000..2f852432 --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test4 @@ -0,0 +1,23 @@ +#!/bin/sh +# Sunearqcsv_test4 - for program SUNEARQCSV +# Author: Andre Latour, Nov 2024 +# + echo "---------------------------------------------------------------------------" + echo "--- Using base1 and target3. " + echo "--- Testing 3 dimensions with third dimension NOT used for nearest, " + echo "--- but used for relative-range-limit. Note that nopoint=0 means that " + echo "--- values for k are -9999 when there is no base point within range AND " + echo "--- several other names will retain their values. " + echo "---------------------------------------------------------------------------" +# + sunearqcsv qin=base1.csv dimx=h dimy=e dimr=c typer=-2 minr=-10 maxr=10 \ + tin=target3.csv timx=xxh timy=xxe timr=xxc \ + qout=qout004.csv nopoint=0 oreps=eye,f,g,k +# + sunearqcsv qin=base1.csv dimx=h dimy=e dimr=c typer=-2 minr=-200 maxr=200 \ + tin=target3.csv timx=xxh timy=xxe timr=xxc \ + qout=qout004a.csv nopoint=0 oreps=eye,f,g,k +# + echo " Note, there should be differences, just easier to check this way... " + diff qout004.csv qout004a.csv +# diff --git a/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test6 b/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test6 new file mode 100644 index 00000000..6861041a --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test6 @@ -0,0 +1,27 @@ +#!/bin/sh +# Sunearqcsv_test6 - for program SUNEARQCSV +# Author: Andre Latour, Nov 2024 +# + echo "---------------------------------------------------------------------------" + echo "--- This is a test of varying-tuples with otuples options. " + echo "--- This copies the tuples from the tin file, nearest qin, or neither. " + echo "--- Note the vuma values are added to the output BEFORE the tuples. " + echo "--- Note the 4th setup, which inputs the same file and just deletes tuples." + echo "---------------------------------------------------------------------------" +# + sunearqcsv qin=vary6.csv dimx=cdp \ + tin=varytin6.csv timx=cdp \ + qout=qout006.csv otuples=1 oreps=vuma formxy=%.15g +# + sunearqcsv qin=vary6.csv dimx=cdp \ + tin=varytin6.csv timx=cdp \ + qout=qout006a.csv otuples=2 oreps=vuma formxy=%.15g +# + sunearqcsv qin=vary6.csv dimx=cdp \ + tin=varytin6.csv timx=cdp \ + qout=qout006b.csv otuples=0 oreps=vuma formxy=%.15g +# + sunearqcsv qin=varytin6.csv dimx=cdp \ + tin=varytin6.csv timx=cdp \ + qout=qout006c.csv otuples=0 formxy=%.15g +# diff --git a/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test7 b/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test7 new file mode 100644 index 00000000..5831cea3 --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test7 @@ -0,0 +1,27 @@ +#!/bin/sh +# Sunearqcsv_test7 - for program SUNEARQCSV +# Author: Andre Latour, Nov 2024 +# + echo "---------------------------------------------------------------------------" + echo "--- This is a test of fixed-tuples with otuples options. " + echo "--- This copies the tuples from the tin file, nearest qin, or neither. " + echo "--- Note that cdpt is replaced. And xxigi,xxigc are added before tuples. " + echo "--- Note the 4th setup, which inputs the same file and just deletes tuples." + echo "---------------------------------------------------------------------------" +# + sunearqcsv qin=fixed7.csv dimx=cdp \ + tin=fixedtin7.csv timx=cdp \ + qout=qout007.csv otuples=1 oreps=cdpt,xxigi,xxigc formxy=%.15g +# + sunearqcsv qin=fixed7.csv dimx=cdp \ + tin=fixedtin7.csv timx=cdp \ + qout=qout007a.csv otuples=2 oreps=cdpt,xxigi,xxigc formxy=%.15g +# + sunearqcsv qin=fixed7.csv dimx=cdp \ + tin=fixedtin7.csv timx=cdp \ + qout=qout007b.csv otuples=0 oreps=cdpt,xxigi,xxigc formxy=%.15g +# + sunearqcsv qin=fixedtin7.csv dimx=cdp \ + tin=fixedtin7.csv timx=cdp \ + qout=qout007c.csv otuples=0 formxy=%.15g +# diff --git a/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test8 b/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test8 new file mode 100644 index 00000000..85f57647 --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test8 @@ -0,0 +1,32 @@ +#!/bin/sh +# Sunearqcsv_test8 - for program SUNEARQCSV +# Author: Andre Latour, Nov 2024 +# + echo "---------------------------------------------------------------------------" + echo "--- This tests combinations of varying-tuples, fixed-tuples, and no-tuples." + echo "---------------------------------------------------------------------------" +# + sunearqcsv qin=fixed7.csv dimx=cdp \ + tin=vary6.csv timx=cdp \ + qout=qout008.csv otuples=2 formxy=%.15g +# + sunearqcsv qin=vary6.csv dimx=cdp \ + tin=fixed7.csv timx=cdp \ + qout=qout008a.csv otuples=2 formxy=%.15g +# + sunearqcsv qin=notuple8.csv dimx=cdp \ + tin=vary6.csv timx=cdp \ + qout=qout008b.csv otuples=1 formxy=%.15g +# + sunearqcsv qin=notuple8.csv dimx=cdp \ + tin=fixed7.csv timx=cdp \ + qout=qout008c.csv otuples=1 formxy=%.15g +# + sunearqcsv qin=vary6.csv dimx=cdp \ + tin=notuple8.csv timx=cdp \ + qout=qout008d.csv otuples=2 formxy=%.15g +# + sunearqcsv qin=fixed7.csv dimx=cdp \ + tin=notuple8.csv timx=cdp \ + qout=qout008e.csv otuples=2 formxy=%.15g +# diff --git a/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test9 b/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test9 new file mode 100644 index 00000000..31d7f1b9 --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/Sunearqcsv_test9 @@ -0,0 +1,29 @@ +#!/bin/sh +# Sunearqcsv_test9 - for program SUNEARQCSV +# Author: Andre Latour, Nov 2024 +# + echo "---------------------------------------------------------------------------" + echo "--- This test uses q-files of velocity stuff from the Sunmocsv tests. " + echo "--- It removes the varying tuples and then puts them back on again. " + echo "--- It removes the fixed tuples and then puts them back on again. " + echo "--- The resulting q-files should have the same values as the originals. " + echo "--- (You can replace qx1 with qoutqx1a and qx5 with qoutqx5a in the " + echo "--- tests within Sunmocsv to confirm they work the same). " + echo "---------------------------------------------------------------------------" +# + sunearqcsv qin=qx1.csv dimx=cdp \ + tin=qx1.csv timx=cdp \ + qout=qoutqx1.csv otuples=0 formxy=%.15g +# + sunearqcsv qin=qx5.csv dimx=cdp \ + tin=qx5.csv timx=cdp \ + qout=qoutqx5.csv otuples=0 formxy=%.15g +# + sunearqcsv qin=qx1.csv dimx=cdp \ + tin=qoutqx1.csv timx=cdp \ + qout=qoutqx1a.csv otuples=2 formxy=%.15g +# + sunearqcsv qin=qx5.csv dimx=cdp \ + tin=qoutqx5.csv timx=cdp \ + qout=qoutqx5a.csv otuples=2 formxy=%.15g +# diff --git a/src/demos/Geom3D/Sunearqcsv/base1.csv b/src/demos/Geom3D/Sunearqcsv/base1.csv new file mode 100644 index 00000000..f1f925a7 --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/base1.csv @@ -0,0 +1,11 @@ +C_SU_SETID,Q,,,,,,,,,,,, +C_SU_NAMES,,,,,,,,,,,,, +C_SU_ID,b,c,d,e,f,g,h,eye,j,k,el,m,n +Q,10,110,210,310,410,510,610,710,880,910,1010,1110,1210 +Q,20,120,220,320,420,520,620,720,870,920,1020,1120,1220 +Q,30,130,230,330,430,530,630,730,860,930,1030,1130,1230 +Q,40,140,240,340,440,540,640,740,850,940,1040,1140,1240 +Q,50,150,250,350,450,550,650,750,840,950,1050,1150,1250 +Q,60,160,260,360,460,560,660,760,830,960,1060,1160,1260 +Q,70,170,270,370,470,570,670,770,820,970,1070,1170,1270 +Q,80,180,280,380,480,580,680,780,810,980,1080,1180,1280 diff --git a/src/demos/Geom3D/Sunearqcsv/base2.csv b/src/demos/Geom3D/Sunearqcsv/base2.csv new file mode 100644 index 00000000..baead032 --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/base2.csv @@ -0,0 +1,11 @@ +C_SU_SETID,Q,,,,,,,,,,,, +C_SU_NAMES,,,,,,,,,,,,, +C_SU_ID,b,c,h,g,eye,el,m,n,d,e,f,j,k +Q,10,110,610,510,710,1010,1110,1210,210,310,410,880,910 +Q,20,120,620,520,720,1020,1120,1220,220,320,420,870,920 +Q,60,160,660,560,760,1060,1160,1260,260,360,460,830,960 +Q,70,170,670,570,770,1070,1170,1270,270,370,470,820,970 +Q,80,180,680,580,780,1080,1180,1280,280,380,480,810,980 +Q,30,130,630,530,730,1030,1130,1230,230,330,430,860,930 +Q,40,140,640,540,740,1040,1140,1240,240,340,440,850,940 +Q,50,150,650,550,750,1050,1150,1250,250,350,450,840,950 diff --git a/src/demos/Geom3D/Sunearqcsv/fixed7.csv b/src/demos/Geom3D/Sunearqcsv/fixed7.csv new file mode 100644 index 00000000..1fa3467f --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/fixed7.csv @@ -0,0 +1,18 @@ +C_SU_SETID,Q +C_SU_FORMS +C_SU_ID,%.20g,%.20g,%.20g,%.20g,%.2f,%.2f,%.2f,%.2f,%.2f,%.2f,%.2f,%.2f +C_SU_NDIMS,tims,,,,,,,,0,190,850,1800 +C_SU_NAMES +C_SU_ID,cdp,cdpt,xxigi,xxigc,gx,gy,sx,sy,vnmo,vnmo,vnmo,vnmo +Q,21,7461,21,1,500.00,0.00,339049.86,5540266.91,1021.00,2000.00,4777.00,5000.00 +Q,41,7921,41,1,1000.00,0.00,339299.72,5539833.81,1041.00,3500.00,4077.00,5000.00 +Q,81,71841,81,1,2000.00,0.00,339799.44,5538967.63,1081.00,3500.00,4000.00,5077.00 +Q,101,72301,101,1,2500.00,0.00,340049.31,5538534.54,1101.00,5500.00,6000.00,8777.00 +Q,1473,7473,21,13,500.00,600.00,339569.57,5540566.74,2473.00,3333.00,5777.00,6666.00 +Q,1493,7933,41,13,1000.00,600.00,339819.43,5540133.65,2493.00,3577.00,4000.00,5000.00 +Q,1533,71853,81,13,2000.00,600.00,340319.16,5539267.46,2777.00,3500.00,4000.00,5000.00 +Q,1553,72313,101,13,2500.00,600.00,340569.02,5538834.37,2553.00,3577.00,4000.00,5000.00 +Q,2683,7483,21,23,500.00,1100.00,340002.67,5540816.60,3683.00,3700.00,4777.00,5000.00 +Q,2703,7943,41,23,1000.00,1100.00,340252.53,5540383.51,3703.00,3800.00,4770.00,5000.00 +Q,2743,71863,81,23,2000.00,1100.00,340752.25,5539517.32,3743.00,3800.00,4070.00,5500.00 +Q,2763,72323,101,23,2500.00,1100.00,341002.11,5539084.23,3763.00,5000.00,6777.00,8000.00 diff --git a/src/demos/Geom3D/Sunearqcsv/fixedtin7.csv b/src/demos/Geom3D/Sunearqcsv/fixedtin7.csv new file mode 100644 index 00000000..231af63b --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/fixedtin7.csv @@ -0,0 +1,16 @@ +C_SU_SETID,Q +C_SU_FORMS +C_SU_ID,%.20g,%.20g,%.20g,%.20g,%.2f,%.2f,%.2f,%.2f,%.2f,%.2f,%.2f,%.2f +C_SU_NDIMS,tnmo,,,,,,,,0,222,888,1800 +C_SU_NAMES +C_SU_ID,cdp,cdpt,igi,igc,gx,gy,sx,sy,vels,vels,vels,vels +Q,21,461,21,1,500.00,0.00,339049.86,5540266.91,1021.00,2000.00,4000.00,5000.00 +Q,41,921,41,1,1000.00,0.00,339299.72,5539833.81,1041.00,3500.00,4000.00,5000.00 +Q,101,2301,101,1,2500.00,0.00,340049.31,5538534.54,1101.00,5500.00,6000.00,8000.00 +Q,1473,473,21,13,500.00,600.00,339569.57,5540566.74,2473.00,3333.00,5555.00,6666.00 +Q,1493,933,41,13,1000.00,600.00,339819.43,5540133.65,2493.00,3500.00,4000.00,5000.00 +Q,1533,1853,81,13,2000.00,600.00,340319.16,5539267.46,2533.00,3500.00,4000.00,5000.00 +Q,2683,483,21,23,500.00,1100.00,340002.67,5540816.60,3683.00,3700.00,4000.00,5000.00 +Q,2703,943,41,23,1000.00,1100.00,340252.53,5540383.51,3703.00,3800.00,4000.00,5000.00 +Q,2743,1863,81,23,2000.00,1100.00,340752.25,5539517.32,3743.00,3800.00,4000.00,5500.00 +Q,2763,2323,101,23,2500.00,1100.00,341002.11,5539084.23,3763.00,5000.00,6000.00,8000.00 diff --git a/src/demos/Geom3D/Sunearqcsv/notuple8.csv b/src/demos/Geom3D/Sunearqcsv/notuple8.csv new file mode 100644 index 00000000..46ee3ef7 --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/notuple8.csv @@ -0,0 +1,15 @@ +C_SU_SETID,Q +C_SU_FORMS +C_SU_ID,%.15g,%.15g,%.15g,%.15g,%.15g,%.15g,%.15g,%.15g +C_SU_NAMES +C_SU_ID,cdp,cdpt,igi,igc,gx,gy,sx,sy +Q,21,461,21,1,500,0,339049.86,5540266.91 +Q,41,921,41,1,1000,0,339299.72,5539833.81 +Q,101,2301,101,1,2500,0,340049.31,5538534.54 +Q,1473,473,21,13,500,600,339569.57,5540566.74 +Q,1493,933,41,13,1000,600,339819.43,5540133.65 +Q,1533,1853,81,13,2000,600,340319.16,5539267.46 +Q,2683,483,21,23,500,1100,340002.67,5540816.6 +Q,2703,943,41,23,1000,1100,340252.53,5540383.51 +Q,2743,1863,81,23,2000,1100,340752.25,5539517.32 +Q,2763,2323,101,23,2500,1100,341002.11,5539084.23 diff --git a/src/demos/Geom3D/Sunearqcsv/qx1.csv b/src/demos/Geom3D/Sunearqcsv/qx1.csv new file mode 100644 index 00000000..96c4a861 --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/qx1.csv @@ -0,0 +1,17 @@ +C_SU_SETID,Q +C_SU_FORMS +C_SU_ID,%.20g,%.20g,%g,%g,%g,%g,%g,%g,%g,%g +C_SU_NAMES +C_SU_ID,cdp,tnmo,vnmo,tnmo,vnmo,tnmo,vnmo,tnmo,vnmo +Q,21,0,1200,190,2000,450,4000,800,5000 +Q,41,0,1800,190,3500,450,4000,800,5000 +Q,81,0,1800,190,3500,450,4000,800,5000 +Q,101,0,2800,190,5500,450,6000,800,8000 +Q,1473,0,3333,200,5555,500,6666 +Q,1493,0,1800,190,3500,450,4000,800,5000 +Q,1533,0,1800,190,3500,450,4000,800,5000 +Q,1553,0,1800,190,3500,450,4000,800,5000 +Q,2683,0,1200,190,2000,450,4000,800,5000 +Q,2703,0,1800,190,3500,450,4000,800,5000 +Q,2743,0,1800,190,3500,450,4000,800,5000,900,5500 +Q,2763,0,2800,190,5500,450,6000,800,8000 diff --git a/src/demos/Geom3D/Sunearqcsv/qx5.csv b/src/demos/Geom3D/Sunearqcsv/qx5.csv new file mode 100644 index 00000000..7b88f8e0 --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/qx5.csv @@ -0,0 +1,21 @@ +C_SU_SETID,Q +C_SU_FORMS +C_SU_ID,%.20g,%.20g,%8.0f,%8.0f,%8.0f,%8.0f +C_SU_NDIMS,tnmo,,0,0.004,0.008,0.012 +C_SU_NAMES +C_SU_ID,cdp,cdpt,vnmo,vnmo,vnmo,vnmo +Q,136,1, 1000136, 2000015, 3000002, 4000000 +Q,172,1, 1000172, 2000051, 3000002, 4000000 +Q,182,1, 1000182, 2000061, 3000002, 4000000 +Q,212,1, 1000212, 2000091, 3000002, 4000000 +Q,232,1, 1000232, 2000111, 3000002, 4000000 +Q,257,1, 1000257, 2000015, 3000003, 4000000 +Q,293,1, 1000293, 2000051, 3000003, 4000000 +Q,303,1, 1000303, 2000061, 3000003, 4000000 +Q,333,1, 1000333, 2000091, 3000003, 4000000 +Q,353,1, 1000353, 2000111, 3000003, 4000000 +Q,620,1, 1000620, 2000015, 3000006, 4000000 +Q,656,1, 1000656, 2000051, 3000006, 4000000 +Q,666,1, 1000666, 2000061, 3000006, 4000000 +Q,696,1, 1000696, 2000091, 3000006, 4000000 +Q,716,1, 1000716, 2000111, 3000006, 4000000 diff --git a/src/demos/Geom3D/Sunearqcsv/target1.csv b/src/demos/Geom3D/Sunearqcsv/target1.csv new file mode 100644 index 00000000..2a8d7512 --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/target1.csv @@ -0,0 +1,11 @@ +C_SU_SETID,Q,,,,,,,,,,,, +C_SU_NAMES,,,,,,,,,,,,, +C_SU_ID,b,c,d,e,f,g,h,eye,j,k,el,m,n +Q,11,111,211,311,411,515,611,711,881,911,1011,1111,1211 +Q,21,121,221,321,421,525,621,721,871,921,1021,1121,1221 +Q,31,131,231,331,431,535,631,731,861,931,1031,1131,1231 +Q,41,141,241,341,441,545,641,741,851,941,1041,1141,1241 +Q,51,151,251,351,451,555,651,751,841,951,1051,1151,1251 +Q,61,161,261,361,461,565,661,761,831,961,1061,1161,1261 +Q,71,171,271,371,471,575,671,771,821,971,1071,1171,1271 +Q,81,181,281,381,481,585,681,781,811,981,1081,1181,1281 diff --git a/src/demos/Geom3D/Sunearqcsv/target2.csv b/src/demos/Geom3D/Sunearqcsv/target2.csv new file mode 100644 index 00000000..b0596064 --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/target2.csv @@ -0,0 +1,9 @@ +C_SU_SETID,Q,,,,,,,,,, +C_SU_NAMES,,,,,,,,,,, +C_SU_ID,xxc,d,b,xxe,f,g,eye,m,n,h,el +Q,111,211,11,311,411,515,711,1111,1211,611,1011 +Q,131,231,31,331,431,535,731,1131,1231,631,1031 +Q,161,261,61,361,461,565,761,1161,1261,661,1061 +Q,171,271,71,371,471,575,771,1171,1271,671,1071 +Q,121,221,21,321,421,525,721,1121,1221,621,1021 +Q,151,251,51,351,451,555,751,1151,1251,651,1051 diff --git a/src/demos/Geom3D/Sunearqcsv/target3.csv b/src/demos/Geom3D/Sunearqcsv/target3.csv new file mode 100644 index 00000000..a5f4dd54 --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/target3.csv @@ -0,0 +1,9 @@ +C_SU_SETID,Q,,,,,,,,,, +C_SU_NAMES,,,,,,,,,,, +C_SU_ID,xxc,d,b,xxe,f,g,eye,m,n,xxh,el +Q,111,211,11,311,411,515,711,1111,1211,611,1011 +Q,231,231,31,331,431,535,731,1131,1231,631,1031 +Q,161,261,61,361,461,565,761,1161,1261,661,1061 +Q,371,271,71,371,471,575,771,1171,1271,671,1071 +Q,121,221,21,321,421,525,721,1121,1221,621,1021 +Q,151,251,51,351,451,555,751,1151,1251,651,1051 diff --git a/src/demos/Geom3D/Sunearqcsv/vary6.csv b/src/demos/Geom3D/Sunearqcsv/vary6.csv new file mode 100644 index 00000000..7255f19b --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/vary6.csv @@ -0,0 +1,17 @@ +C_SU_NAMES +C_SU_ID,cdp,cdpt,numa,vuma,tims,tupb,vels,tims,tupb,vels +Q,136,1,3136.000,1136.000,0.000,15136.000,10136.000,0.004,25015.000,20015.000,0.008,35002.000,30002.000,0.012,45006.666,40000.000 +Q,172,2,3172.000,1172.000,0.000,15172.000,10172.000,0.004,25051.000,20051.000,0.008,35002.000,30006.666 +Q,182,3,3182.000,1182.000,0.000,15182.000,10182.000,0.004,25061.000,20061.000,0.008,35002.000,30002.000,0.012,45006.666,40000.000 +Q,212,4,3212.000,1212.000,0.000,15212.000,10212.000,0.004,25091.000,20091.666 +Q,232,5,3232.000,1232.000,0.000,15232.000,10232.000,0.004,25111.000,20111.666,,, +Q,257,6,3257.000,1257.000,0.000,15257.000,10257.000,0.004,25015.000,20015.000,0.008,35003.000,30003.666,0.012,45006.000,40000.000 +Q,293,7,3293.000,1293.000,0.000,15293.000,10293.000,0.004,25051.000,20051.000,0.008,35003.000,30003.666,0.012,45006.000,40000.000 +Q,303,8,3303.000,1303.000,0.000,15303.000,10303.000,0.004,25061.000,20061.000,0.008,35003.000,30003.666,0.012,45006.000,40000.000,0.2,55000,50000 +Q,333,9,3333.000,1333.000,0.000,15333.000,10333.000,0.004,25091.000,20091.000,0.008,35003.000,30003.666,0.012,45006.000,40000.000 +Q,353,10,3353.000,1353.000,0.000,15353.000,10353.000,0.004,25111.000,20111.000,0.008,35003.000,30003.666,0.012,45006.000,40000.000 +Q,620,11,3620.000,1620.000,0.000,15620.000,10620.000,0.004,25015.000,20015.000,0.008,35006.000,30006.666,0.012,45006.000,40000.000,,, +Q,656,12,3656.000,1656.000,0.000,15656.000,10656.000,0.004,25051.000,20051.000,0.008,35006.000,36666.000,0.012,45006.000,40000.000 +Q,666,13,3666.000,1666.000,0.000,15666.000,10666.000,0.004,25061.000,20061.000,0.008,35006.000,36666.000,0.012,45006.000,40000.000,,,,,, +Q,696,14,3696.000,1696.000,0.000,15696.000,10696.000,0.004,25091.000,20091.000,0.008,35006.000,30606.000,0.012,45006.000,40000.000 +Q,716,15,3716.000,1716.000,0.000,15716.000,10716.000,0.004,25111.000,20111.000,0.008,35006.000,30066.000,0.012,45006.000,40000.000 diff --git a/src/demos/Geom3D/Sunearqcsv/varytin6.csv b/src/demos/Geom3D/Sunearqcsv/varytin6.csv new file mode 100644 index 00000000..d77216af --- /dev/null +++ b/src/demos/Geom3D/Sunearqcsv/varytin6.csv @@ -0,0 +1,15 @@ +C_SU_NAMES +C_SU_ID,cdp,numa,xvuma,tnmo,dpth,vnmo,tnmo,dpth,vnmo +Q,136,3136.000,6136.000,0.000,15136.000,10136.000,0.004,25015.000,20015.000,0.008,35002.000,30002.000,0.012,45000.000,40000.000 +Q,172,3172.000,6172.000,0.000,15172.000,10172.000,0.004,25051.000,20051.000,0.008,35002.000,30002.000 +Q,182,3182.000,6182.000,0.000,15182.000,10182.000,0.004,25061.000,20061.000,0.008,35002.000,30002.000,0.012,45000.000,40000.000 +Q,212,3200.000,6200.000,0.000,15000.000,10000.000,0.004,25000.000,20000.000 +Q,232,3232.000,6232.000,0.000,15232.000,10232.000,0.004,25111.000,20111.000,,, +Q,257,3257.000,6257.000,0.000,15257.000,10257.000,0.004,25015.000,20015.000,0.008,35003.000,30003.000,0.012,45000.000,40000.000 +Q,293,3293.000,6293.000,0.000,15293.000,10293.000,0.004,25051.000,20051.000,0.008,35003.000,30003.000,0.012,45000.000,40000.000 +Q,333,3333.000,6333.000,0.000,15333.000,10333.000,0.004,25091.000,20091.000,0.008,35003.000,30003.000,0.012,45000.000,40000.000 +Q,353,3353.000,6353.000,0.000,15353.000,10353.000,0.004,25111.000,20111.000,0.008,35003.000,30003.000,0.012,45000.000,40000.000 +Q,656,3656.000,6656.000,0.000,15656.000,10656.000,0.004,25051.000,20051.000,0.008,35006.000,30006.000,0.012,45000.000,40000.000 +Q,666,3666.000,6666.000,0.000,15666.000,10666.000,0.004,25061.000,20061.000,0.008,35006.000,30006.000,0.012,45000.000,40000.000,,,,,, +Q,696,3696.000,6696.000,0.000,15696.000,10696.000,0.004,25091.000,20091.000,0.008,35006.000,30006.000,0.012,45000.000,40000.000 +Q,716,3716.000,6716.000,0.000,15716.000,10716.000,0.004,25111.000,20111.000,0.008,35006.000,30006.000,0.012,45000.000,40000.000 From 4fa3035c7d4e357955776997c810a6907224b84e Mon Sep 17 00:00:00 2001 From: l2r <80277078+l2r-git@users.noreply.github.com> Date: Sat, 9 Nov 2024 11:28:42 -0800 Subject: [PATCH 5/5] New program sunearqcsv. sunearqcsv copies nearest values from one q-file to another. --- src/su/main/headers/Makefile | 1 + src/su/main/headers/sunearqcsv.c | 1518 ++++++++++++++++++++++++++++++ 2 files changed, 1519 insertions(+) create mode 100644 src/su/main/headers/sunearqcsv.c diff --git a/src/su/main/headers/Makefile b/src/su/main/headers/Makefile index 8ff2c02f..ab25d3fb 100644 --- a/src/su/main/headers/Makefile +++ b/src/su/main/headers/Makefile @@ -36,6 +36,7 @@ PROGS = \ $B/sulcthw \ $B/sulhead \ $B/sunearcsv \ + $B/sunearqcsv \ $B/sunomstreamers \ $B/supaste \ $B/suprofcsv \ diff --git a/src/su/main/headers/sunearqcsv.c b/src/su/main/headers/sunearqcsv.c new file mode 100644 index 00000000..28f447fb --- /dev/null +++ b/src/su/main/headers/sunearqcsv.c @@ -0,0 +1,1518 @@ +/* Copyright (c) Colorado School of Mines, 2022.*/ +/* All rights reserved. */ + +/* SUNEARCSV: $Revision: 1.01 $ ; $Date: 2024/10/20 00:00:01 $ */ + +#include "su.h" +#include "segy.h" +#include +#include "qdefine.h" + +/*********************** self documentation ******************************/ +char *sdoc[] = { +" ", +" SUNEARQCSV - Assign Values To Target Q-file From Nearest In Another Q-file.", +" ", +" sunearqcsv [parameters]. (No traces in or out). ", +" ", +" Parameters: ", +" ", +" qin= Input q-file to copy nearest values FROM. The values to copy are ", +" specified via the oreps parameter and the otuples parameter. ", +" ", +" dimx= Dimension name X. Name must exist in the qin q-file. ", +" dimy= Dimension name Y. Name must exist in the qin q-file. ", +" dimr= Dimension name R. Name must exist in the qin q-file. ", +" dima= Dimension name A. Name must exist in the qin q-file. ", +" dimb= Dimension name B. Name must exist in the qin q-file. ", +" dimc= Dimension name C. Name must exist in the qin q-file. ", +" dimd= Dimension name D. Name must exist in the qin q-file. ", +" dime= Dimension name E. Name must exist in the qin q-file. ", +" dimf= Dimension name F. Name must exist in the qin q-file. ", +" Note: No defaults. At least 1 dimension must be used. And for every ", +" qin dimension parameter name specified the corresponding ", +" tin target dimension parameter name must be specified. ", +" ", +" The following 3 parameters can be used for any of the dimensions. ", +" Just substitute x,y,r,a,b,c,d,e,f as the ending character, ", +" for instance: typer,minr,maxr for dimension R. ", +" Typically these 3 are used just to specify an additional dimension that is ", +" used to restrict the search range (see typer options -1 and -2 below). ", +" ", +" typer=1 Use in Pythagorean Nearest (squared difference to tin location), ", +" and minr,maxr are unchanging extent ranges for this dimension. ", +" =2 Use in Pythagorean Nearest and minr,maxr are the relative extent ", +" range to include for this dimension (the input tin value for ", +" this dimension is added to minr,maxr to form the extent range). ", +" =-1 Do not use in Pythagorean Nearest, and minr,maxr are unchanging ", +" extent ranges for this dimension. ", +" =-2 Do not use in Pythagorean Nearest and minr,maxr are the relative ", +" extent range to include for this dimension (the input tin value ", +" for this dimension is added to minr,maxr to form extent range). ", +" Note: Negative types means the difference between the points and ", +" the tin is NOT added to Pythagorean sum for specified dimension(s). ", +" So, the nearest point is determined as if that dimension was NOT ", +" specified at all. However, the min,max for that dimension are still ", +" used. This results in finding the nearest point considering only the ", +" type>0 dimensions, but restricted by ranges of type<0 dimension(s). ", +" Type=-2 can be used for situations where a profile approaches itself ", +" or intersects itself or overlaps itself. In those case, a third ", +" dimension value (such as station number) can be used to restrict the ", +" search range for each tin to only the part of the profile with ", +" approximately the same station numbers as the tin midpoint station. ", +" minr= Extent Range Min. Typically negative. Default is no min limit. ", +" Greater or equal to this value is in range. ", +" maxr= Extent Range Max. Typically positive. Default is no max limit. ", +" Strictly less than this value is in range. The min must be less ", +" than max, but the range does not have to be symmetric or centred ", +" (that is, ranges such as min=-1000 and max=-200 are allowed). ", +" ", +" ordp=asis Point Order Name. Sometimes there are two-or-more qin points ", +" which are equally near the tin target location. In those cases, ", +" the nearest point is set to the higher ordered point. ", +" The default (asis) means use the input order of qin records. ", +" Specifying a name here will cause the input qin records to be ", +" sorted internally by the values of the specified name. ", +" Note: The sunearcsv program has a similar parameter called keyp. ", +" This parameter matches its ordering functionality. ", +" ", +" tin= Target Input q-file to copy nearest values FROM qin q-file. ", +" ", +" timx= Target Dimension name X. Name must exist in the tin q-file. ", +" timy= Target Dimension name Y. Name must exist in the tin q-file. ", +" timr= Target Dimension name R. Name must exist in the tin q-file. ", +" tima= Target Dimension name A. Name must exist in the tin q-file. ", +" timb= Target Dimension name B. Name must exist in the tin q-file. ", +" timc= Target Dimension name C. Name must exist in the tin q-file. ", +" timd= Target Dimension name D. Name must exist in the tin q-file. ", +" time= Target Dimension name E. Name must exist in the tin q-file. ", +" timf= Target Dimension name F. Name must exist in the tin q-file. ", +" Note: A target dimension must be specified for each dim dimension used.", +" ", +" qout=qout.csv Output q-file name. ", +" ", +" oreps= Output Replace Name List. ", +" Names listed here must be in the non-tuples part of qin q-file. ", +" If name here already exists in non-tuples part of tin file, its ", +" nearest value from the qin file are REPLACED in the output file. ", +" If name does not exist in non-tuples part of tin file, its name ", +" and nearest value from the qin file are ADDED to the output file.", +" Names and values NOT listed here are COPIED from tin to output. ", +" Note: This parameter does not have to be specified but the resulting ", +" output might just be a copy of input depending on the following ", +" otuples parameter and whether tuples exist in qin or tin. ", +" ", +" otuples=1 Output the tuples asis from tin file (if they exist). ", +" =0 Do not output any tuples (even if they exist). ", +" =2 Output the tuples of nearest point in qin file (if they exist). ", +" Tuples cannot be intermixed (either all from tin, or from qin). ", +" Tuples can have varying or fixed amounts (see Q File Standards). ", +" Tuple names copied from the qin file are NOT error-checked to see", +" if the same name exists in the non-tuple names in the tin file. ", +" ", +" nopoint=1 If no qin point is found within the extents, error-halt (before ", +" outputting the first q-record that could not find a qin point). ", +" =0 Continue without any update of oreps values for q-records that ", +" could not find a qin point in extents (print their count at end).", +" *** Since the default extent min,max are unlimited, this option has ", +" no effect unless you specify them. ", +" Note: If you choose to use option 0, then values for oreps that exist ", +" in tin file are output with their tin values. But if the oreps ", +" name does not exist in the tin file, then output value is -9999. ", +" ", +" formxy=%.20g The C format code for printing all values to qout records. ", +" Note that the default format prints up to 20 digits ", +" (but not trailing zeroes to the right of the decimal point). ", +" ", +" ", +" The following 3 parameters affect cpu time, but not results. The search ", +" is done by building a kdtree from the dimension values in the qin file. ", +" The code that builds the kdtree is reasonably standard and simplistic. ", +" But the kdtree search code is slightly unusual due to some typer options. ", +" ", +" sfunc=2 Search via the cycle_for_near function. This option is usually ", +" fastest. This option uses the sdist and smult parameters. ", +" =1 Search via the find_near function. This option may be faster ", +" if you are specifying small extent ranges. ", +" This option does not use the sdist and smult parameters. ", +" =0 Search via the brute_near function. This option may be faster ", +" if there are only a small number of points in the qin file. ", +" This option does not use the sdist and smult parameters. ", +" (To keep the code simpler, this option still allocates and ", +" builds the kdtree. But the tree is not used for searching). ", +" =-1 Search via the find_near function and confirm the results using ", +" the brute_near function. This tests the find_near function and ", +" can also be used to determine the speed difference between them. ", +" This option does not use the sdist and smult parameters. ", +" =-2 Search via the cycle_for_near function and confirm the results ", +" using brute_near function. This tests cycle_for_near function ", +" and can also be used to determine speed difference between them. ", +" This option uses the sdist and smult parameters. ", +" sdist=100 Initial search distance. This parameter is only used ", +" if sfunc=2 or -2. If positive, this value is added to the ", +" distance between the previous tin and its nearest qin point ", +" and used as the initial search distance for the current tin. ", +" If negative, the initial search distance for all tin is set ", +" to this absolute value. ", +" smult=2 Search multiplier. This parameter is only used if sfunc=2 or -2. ", +" If the nearest qin point is not found after searching with the ", +" initial search distance, the search distance is multiplied by ", +" this value, and search is performed again. This repeats until ", +" finding the nearest point (or all min,max ranges are exceeded). ", +" ", +" ------------------------------------------------------------------ ", +" ------------------------------------------------------------------ ", +" ", +NULL}; + +/* Created: Oct 2024: Andre Latour */ +/* This program started from sunearcsv. */ +/**************** end self doc *******************************************/ + +segy tr; + +struct QInfo *RecInfo; /* Storage for all function location value pointers */ +struct QInfo *RecInfot; /* Storage for all function location value pointers */ +int locp = -1; + +int compSort1 (const void * q1, const void * q2) ; /* comparison function for qsort */ + +/* Note: I make no claim that this is a particularly good kd tree implementation. */ +/* It is not explicitly balanced (it has an option to get approximate balance). */ +/* Note: Option tree_nt is unlikely to exist in other kd tree implementations. */ +/* It exists because some crooked-profiles (land) or coil-profiles (marine) */ +/* curve back-over-top-of-themselves. These self-intersections mean the profile */ +/* is not a function (in the mathematical sense). That is, just considering XYs */ +/* a tin midpoint can get confused as to which part of the profile it should */ +/* belong to. Similar issues arise if the profile just curves back NEAR itself. */ + +typedef struct node { + unsigned long long elem; + struct node * l; + struct node * r; +} node; + +/* Note: Yes, I could also have made a structure named tree, and then passed it into */ +/* the functions instead of some of the individual arguments. That would have */ +/* reduced the function arguments, but make it more confusing for new coders. */ + +void connect_nodes (node *tree_nodes, unsigned long long tree_numc, double **tree_dl, int tree_numd, + int ihop); + +void connect_nodes_your (node *tree_nodes, unsigned long long tree_numc, double **tree_dl, int tree_numd, + unsigned long long *your_order); + +void connect_all (node *tree_nodes, unsigned long long tree_numc, double **tree_dl, int tree_numd); + +void find_in (node *tree_nodes, double **tree_dl, int tree_numd, + double *extent_min, double *extent_max, + unsigned long long *out_elem, unsigned long long *num_out); + +void find_in_rcur (double **tree_dl, int tree_numd, node *now_node, int naxe, + double *extent_min, double *extent_max, + unsigned long long *out_elem, unsigned long long *num_out); + +void find_near (node *tree_nodes, double **tree_dl, int *tree_nt, int tree_numd, + double *extent_min, double *extent_max, double *target, + unsigned long long *near_elem, double *near_dist, unsigned long long *num_found); + +void find_near_rcur (double **tree_dl, int *tree_nt, int tree_numd, node *now_node, int naxe, + double *extent_min, double *extent_max, double *target, + unsigned long long *near_elem, double *near_dist, unsigned long long *num_found); + +void cycle_for_near (node *tree_nodes, double **tree_dl, int *tree_nt, int tree_numd, + double *extent_min, double *extent_max, double *target, + double init_rad, double rad_scal, + unsigned long long *near_elem, double *near_dist, + unsigned long long *num_found, int *ncycles); + +void brute_near (double **tree_dl, unsigned long long tree_numc, int *tree_nt, int tree_numd, + double *extent_min, double *extent_max, double *target, + unsigned long long *near_elem, double *near_dist, unsigned long long *num_found); + +/*----------------------------------------------------------------------*/ + +int main(int argc, char **argv) { + + int ncdp = 0; + int ifixd = 0; /* flag for all tuples same size or vary */ + int iztuple = 0; /* element number where first tuple exists */ + int ktuple = 0; /* type of tuples (2=pairs, 3=triplets) */ + cwp_String Pname=NULL; /* text file name for Q input file */ + FILE *fpP=NULL; /* file pointer for Q input file */ + + cwp_String *pname = NULL; + cwp_String *ndims = NULL; + int numpname = 0; + double *pindepa = NULL; + int numdind = 0; + + + int ncdpt = 0; + int ifixdt = 0; /* flag for all tuples same size or vary */ + int iztuplet = 0; /* element number where first tuple exists */ + int ktuplet = 0; /* type of tuples (2=pairs, 3=triplets) */ + cwp_String Tname=NULL; /* text file name for Q input file */ + FILE *fpT=NULL; /* file pointer for Q input file */ + + + cwp_String *pnamet = NULL; + cwp_String *ndimst = NULL; + int numpnamet = 0; + double *pindepat = NULL; + int numdindt = 0; + int moret = 0; + + cwp_String formxyt=NULL; + cwp_String formxy=NULL; + cwp_String formxylong=NULL; + int lenformxy = 0; + + int jcdp = 0; + int i = 0; + int j = 0; + int k = 0; + int errwarn = 0; + + cwp_String keyn[9]; + int locn[9]; + int tree_numd = 0; + + cwp_String keyt[9]; + int locnt[9]; + + int num_ordp = 0; + + int num_oreps = 0; + cwp_String *oreps = NULL; + int *oloc = NULL; + int *oloct = NULL; + int *olocm = NULL; + + int otuples = 1; + int nopoint = 1; + double dmiss = -9999.; + + double sdist = 100.0; + double smult = 2.; + int ncycles = 0; + int tcycles = 0; + int sfunc = 2; + int scheck = 0; + int ncheck = 0; + int sdadd = 1; + double sdist2 = 100.0; + + unsigned long long near_elem = 0; + unsigned long long num_found = 0; + double near_dist = 0.; + unsigned long long near_elemc = 0; + unsigned long long num_foundc = 0; + double near_distc = 0.; + int ihop = 1; + int nproct = 0; + int noupdate = 0; + + int *in_rel_extent = NULL; + double *extent_in_min = NULL; + double *extent_in_max = NULL; + double **tree_dl = NULL; + int *tree_nt = NULL; + double *extent_min = NULL; + double *extent_max = NULL; + double *target = NULL; + unsigned long long tree_ncdp = 0; + node *tree_nodes = NULL; + + cwp_String Oname=NULL; /* text file name for O output file */ + FILE *fpO=NULL; /* file pointer for O output file */ + +/* hook up getpar */ + initargs(argc, argv); + requestdoc(1); + + if(isatty(STDIN_FILENO)!=1 || isatty(STDOUT_FILENO)!=1) + err("**** Error: this program does not input or output traces."); + +/* Maximum number of dimensions is 9. */ + + in_rel_extent = ealloc1int(9); + extent_in_min = ealloc1double(9); + extent_in_max = ealloc1double(9); + tree_nt = ealloc1int(9); + extent_min = ealloc1double(9); + extent_max = ealloc1double(9); + target = ealloc1double(9); + tree_dl = ealloc1(9,sizeof(double *)); + +/* Set defaults for extent ranges. */ + + for(i=0; i<9; i++) { + in_rel_extent[i] = 1; /* use in Pythagorean Nearest, has unchanging range */ + extent_in_min[i] = -DBL_MAX; + extent_in_max[i] = DBL_MAX; + } + +/* Note here that keyn is loaded with key names in the order in which the */ +/* dimension parameters are checked next. And keyn is never re-arranged. */ +/* This means that, eventually, the target array (named target) also must */ +/* have its values in the same order. */ + + if(countparval("dimx") > 0) { + getparstring("dimx",keyn+tree_numd ); + if(countparval("typex")>0) getparint("typex",in_rel_extent+tree_numd); + if(countparval("minx") >0) getpardouble("minx",extent_in_min+tree_numd); + if(countparval("maxx") >0) getpardouble("maxx",extent_in_max+tree_numd); + if(countparval("timx") > 0) getparstring("timx",keyt+tree_numd ); + else err("**** Error: Parameter dimx specified, but not timx."); + tree_numd++; + } + else if(countparval("typex")>0 || countparval("minx")>0 || countparval("maxx")>0 || countparval("timx")>0) { + err("**** Error: Parameter typex, minx, maxx, or timx specified, but not dimx."); + } + + if(countparval("dimy") > 0) { + getparstring("dimy",keyn+tree_numd ); + if(countparval("typey")>0) getparint("typey",in_rel_extent+tree_numd); + if(countparval("miny") >0) getpardouble("miny",extent_in_min+tree_numd); + if(countparval("maxy") >0) getpardouble("maxy",extent_in_max+tree_numd); + if(countparval("timy") > 0) getparstring("timy",keyt+tree_numd ); + else err("**** Error: Parameter dimy specified, but not timy."); + tree_numd++; + } + else if(countparval("typey")>0 || countparval("miny")>0 || countparval("maxy")>0 || countparval("timy")>0) { + err("**** Error: Parameter typey, miny, maxy, timy specified, but not dimy."); + } + + if(countparval("dimr") > 0) { + getparstring("dimr",keyn+tree_numd ); + if(countparval("typer")>0) getparint("typer",in_rel_extent+tree_numd); + if(countparval("minr") >0) getpardouble("minr",extent_in_min+tree_numd); + if(countparval("maxr") >0) getpardouble("maxr",extent_in_max+tree_numd); + if(countparval("timr") > 0) getparstring("timr",keyt+tree_numd ); + else err("**** Error: Parameter dimr specified, but not timr."); + tree_numd++; + } + else if(countparval("typer")>0 || countparval("minr")>0 || countparval("maxr")>0 || countparval("timr")>0) { + err("**** Error: Parameter typer, minr, maxr, timr specified, but not dimr."); + } + + if(countparval("dima") > 0) { + getparstring("dima",keyn+tree_numd ); + if(countparval("typea")>0) getparint("typea",in_rel_extent+tree_numd); + if(countparval("mina") >0) getpardouble("mina",extent_in_min+tree_numd); + if(countparval("maxa") >0) getpardouble("maxa",extent_in_max+tree_numd); + if(countparval("tima") > 0) getparstring("tima",keyt+tree_numd ); + else err("**** Error: Parameter dima specified, but not tima."); + tree_numd++; + } + else if(countparval("typea")>0 || countparval("mina")>0 || countparval("maxa")>0 || countparval("tima")>0) { + err("**** Error: Parameter typea, mina, maxa, or tima specified, but not dima."); + } + + if(countparval("dimb") > 0) { + getparstring("dimb",keyn+tree_numd ); + if(countparval("typeb")>0) getparint("typeb",in_rel_extent+tree_numd); + if(countparval("minb") >0) getpardouble("minb",extent_in_min+tree_numd); + if(countparval("maxb") >0) getpardouble("maxb",extent_in_max+tree_numd); + if(countparval("timb") > 0) getparstring("timb",keyt+tree_numd ); + else err("**** Error: Parameter dimb specified, but not timb."); + tree_numd++; + } + else if(countparval("typeb")>0 || countparval("minb")>0 || countparval("maxb")>0 || countparval("timb")>0) { + err("**** Error: Parameter typeb, minb, maxb, or timb specified, but not dimb."); + } + + if(countparval("dimc") > 0) { + getparstring("dimc",keyn+tree_numd ); + if(countparval("typec")>0) getparint("typec",in_rel_extent+tree_numd); + if(countparval("minc") >0) getpardouble("minc",extent_in_min+tree_numd); + if(countparval("maxc") >0) getpardouble("maxc",extent_in_max+tree_numd); + if(countparval("timc") > 0) getparstring("timc",keyt+tree_numd ); + else err("**** Error: Parameter dimc specified, but not timc."); + tree_numd++; + } + else if(countparval("typec")>0 || countparval("minc")>0 || countparval("maxc")>0 || countparval("timc")>0) { + err("**** Error: Parameter typec, minc, maxc, or timc specified, but not dimc."); + } + + if(countparval("dimd") > 0) { + getparstring("dimd",keyn+tree_numd ); + if(countparval("typed")>0) getparint("typed",in_rel_extent+tree_numd); + if(countparval("mind") >0) getpardouble("mind",extent_in_min+tree_numd); + if(countparval("maxd") >0) getpardouble("maxd",extent_in_max+tree_numd); + if(countparval("timd") > 0) getparstring("timd",keyt+tree_numd ); + else err("**** Error: Parameter dimd specified, but not timd."); + tree_numd++; + } + else if(countparval("typed")>0 || countparval("mind")>0 || countparval("maxd")>0 || countparval("timd")>0) { + err("**** Error: Parameter typed, mind, maxd, timd specified, but not dimd."); + } + + if(countparval("dime") > 0) { + getparstring("dime",keyn+tree_numd ); + if(countparval("typee")>0) getparint("typee",in_rel_extent+tree_numd); + if(countparval("mine") >0) getpardouble("mine",extent_in_min+tree_numd); + if(countparval("maxe") >0) getpardouble("maxe",extent_in_max+tree_numd); + if(countparval("time") > 0) getparstring("time",keyt+tree_numd ); + else err("**** Error: Parameter dime specified, but not time."); + tree_numd++; + } + else if(countparval("typee")>0 || countparval("mine")>0 || countparval("maxe")>0 || countparval("time")>0) { + err("**** Error: Parameter typee, mine, maxe, or time specified, but not dime."); + } + + if(countparval("dimf") > 0) { + getparstring("dimf",keyn+tree_numd ); + if(countparval("typef")>0) getparint("typef",in_rel_extent+tree_numd); + if(countparval("minf") >0) getpardouble("minf",extent_in_min+tree_numd); + if(countparval("maxf") >0) getpardouble("maxf",extent_in_max+tree_numd); + if(countparval("timf") > 0) getparstring("timf",keyt+tree_numd ); + else err("**** Error: Parameter dimf specified, but not timf."); + tree_numd++; + } + else if(countparval("typef")>0 || countparval("minf")>0 || countparval("maxf")>0 || countparval("timf")>0) { + err("**** Error: Parameter typef, minf, maxf, timf specified, but not dimf."); + } + + if(tree_numd<1) err("**** Error: At least 1 dimension name must be specified."); + +/* Resolve a few things. ---------------------------------------------------- */ + + for(i=0; i0) { + tree_nt[i] = 1; /* use as part of Pythagorean Nearest */ + } + else { + tree_nt[i] = 0; /* do not use as part of Pythagorean Nearest */ + in_rel_extent[i] = 0 - in_rel_extent[i]; /* reset to positive 1 or 2 */ + } + + } /* end of for(i=0; i 0) { + getparstring("ordp", &ordp); + if(strcmp(ordp,"asis")!=0) num_ordp = 1; + } + + if(countparval("oreps")>0) { + num_oreps = countparval("oreps"); + oreps = ealloc1(num_oreps,sizeof(cwp_String *)); + getparstringarray("oreps", oreps); + oloc = ealloc1int(num_oreps); + oloct = ealloc1int(num_oreps); + olocm = ealloc1int(num_oreps); + } + + if(!getparint("otuples",&otuples)) otuples = 1; + if(otuples>2 || otuples<0) err("**** Error: otuples must be 0, 1, or 2."); + + if(!getparint("nopoint",&nopoint)) nopoint = 1; + if(nopoint>1 || nopoint<0) err("**** Error: nopoint must be 1 or 0."); + + getparstring("formxy",&formxyt); + if(formxyt==NULL) { + lenformxy = 5; + formxy = ealloc1(lenformxy,1); + strcpy(formxy,"%.20g"); + } + else { + lenformxy = strlen(formxyt); + formxy = ealloc1(lenformxy,1); + strcpy(formxy,formxyt); + } + + formxylong = ealloc1(1+lenformxy,1); + strcpy(formxylong,","); + strcpy(formxylong+1,formxy); + + if(!getparint("sfunc",&sfunc)) sfunc = 2; + if(sfunc<-2 || sfunc>2) err("**** Error: sfunc option is out-of-range."); + + scheck = 0; + if(sfunc<0) { + sfunc = 0 - sfunc; + scheck = 1; + } + + if(!getpardouble("sdist",&sdist)) sdist = 100.; + if(sdist==0.) err("**** Error: sdist cannot be 0."); + + if(!getpardouble("smult",&smult)) smult = 2.; + if(smult<0.) err("**** Error: smult must be non-negative."); + + if(sdist<0.) { + sdist = 0. - sdist; + sdadd = 0; + } + sdist2 = sdist; + +/*-------------------------------------------------------------------------- */ + + getparstring("qin", &Pname); + + fpP = fopen(Pname, "r"); + if(fpP==NULL) err("error: qin input q-file did not open correctly."); + +/* Set input numpname,pname to just store what is going to be accessed. */ +/* numpname>0 is a flag to ONLY store values if they are on pname list. */ +/* numpname<1 is a flag to NOT store values if they are on pname list. */ + +/*otuples=1 Output the tuples asis from tin file. */ +/* =2 Output the tuples of nearest point in qin file. */ +/* =0 Do not output any tuples. */ + + pname = ealloc1(tree_numd + num_ordp + num_oreps,sizeof(cwp_String *)); + numpname = 0; + + if(otuples!=2) { /* Output tuples from qin. So input all values from qin. */ + for (i=0; i100) + err("getviaqfile error: qin file, record %d (wrong comma count, damaged, non-numbers, ...)",errwarn-99); + else if(errwarn>0) err("getviaqfile error: qin file, unrecognized error code %d",errwarn); + + getparstring("tin", &Tname); + + fpT = fopen(Tname, "r"); + if(fpT==NULL) err("error: tin input q-file did not open correctly."); + + errwarn = 0; + getviaqfile(fpT, &pnamet, &numpnamet, &iztuplet, numdindt, + &ktuplet, &ifixdt, &RecInfot, &ncdpt, + &pindepat, &ndimst, &errwarn) ; + + if(errwarn==1) err("getqinfo error: tin file, extra C_SU_NAMES record"); + else if(errwarn==2) err("getqinfo error: tin file, extra C_SU_NDIMS record"); + else if(errwarn==3) err("getqinfo error: tin file, C_SU_ID record not found immediately after C_SU_NAMES record."); + else if(errwarn==11) + err("readqhead error: tin file, if C_SU_NDIMS not vary, its numbers must align with C_SU_NAMES"); + else if(errwarn==12) + err("readqhead error: tin file, C_SU_ID record not found immediately after C_SU_NAMES record."); + else if(errwarn==22) err("getviaqfile error: tin file, C_SU_NDIMS record not same length as C_SU_NAMES record."); + else if(errwarn==23) err("getviaqfile error: tin file, C_SU_NAMES tupled names out-of-order, changed"); + else if(errwarn==24) err("getviaqfile error: tin file, C_SU_NDIMS blank where valid number expected"); + else if(errwarn==25) err("getviaqfile error: tin file, C_SU_NDIMS non-number where valid number expected"); + else if(errwarn==26) err("getviaqfile error: tin file, C_SU_NDIMS value must be same for all members of tuple"); + else if(errwarn==27) err("getviaqfile error: tin file, C_SU_NAMES record followed by C_SU_ID record not found."); + else if(errwarn>100) + err("getviaqfile error: tin file, record %d (wrong comma count, damaged, non-numbers, ...)",errwarn-99); + else if(errwarn>0) err("getviaqfile error: tin file, unrecognized error code %d",errwarn); + + getparstring("qout", &Oname); + + if(strcmp(Pname,Oname) == 0 || strcmp(Tname,Oname) == 0 ) + err("**** Error: qout output file name must be different than qin and tin input file names."); + + fpO = fopen(Oname, "w"); + if (fpO == NULL) err("qfile error: qout output file did not open correctly."); + + checkpars(); + + if(otuples==1 && ifixdt==2) otuples = 0; + if(otuples==2 && ifixd ==2) otuples = 0; + +/*-------------------------------------------------------------------------- */ + + for (i=0; i0) { + locp = -1; + for (i=0; i0) sdist2 = sdist + sqrt(near_dist); + cycle_for_near(tree_nodes,tree_dl,tree_nt,tree_numd, + extent_min, extent_max, target, + sdist2, smult, + &near_elem, &near_dist, &num_found, &ncycles); + tcycles += ncycles; + } + else if(sfunc==1) { + find_near(tree_nodes,tree_dl,tree_nt,tree_numd, + extent_min, extent_max, target, + &near_elem, &near_dist, &num_found); + } + else { + brute_near(tree_dl,tree_ncdp,tree_nt,tree_numd, + extent_min, extent_max, target, + &near_elem, &near_dist, &num_found); + } + + if(scheck==1) { + brute_near(tree_dl,tree_ncdp,tree_nt,tree_numd, + extent_min, extent_max, target, + &near_elemc, &near_distc, &num_foundc); + if(near_elemc!=near_elem || near_distc!=near_dist || num_foundc!=num_found) { + ncheck++; + if(ncheck<10) { + warn("near point: brute=%ld tree=%ld dist diff (squared)=%g equally near: brute=%ld tree=%ld tin counter=%d ", + near_elemc,near_elem,near_distc-near_dist,num_foundc,num_found,nproct+1); + } + } + } + + if(num_found<1) { + if(nopoint>0) err("error: No qin point within extents for tin q-record %d and nopoint option is 1.",nproct+1); + noupdate++; + } + else { + for (i=0; i -1) RecInfot[jcdp].dlots[oloct[i]] = RecInfo[near_elem].dlots[oloc[i]]; + } + } + + nproct++; + + fprintf(fpO,"Q"); + for (i=0; i0) { + for (i=0; i=0; k--) { + fprintf(fpO,formxylong,RecInfot[jcdp].dlots[iztuplet+k*RecInfot[jcdp].nto+i]); + } + } + } + else if(otuples==2) { + if(num_found>0) { + for(i=0; i=0; k--) { + fprintf(fpO,formxylong,RecInfo[near_elem].dlots[iztuple+k*RecInfo[near_elem].nto+i]); + } + } + } + else { + for(k=ktuple-1; k>=0; k--) fprintf(fpO,formxylong,dmiss); + } + } + + fprintf(fpO,"\n"); + + } /* end of for (jcdp=0; jcdp0) warn("Number of q-records not updated due to no qin point within extents = %d (permitted when nopoint=0).",noupdate); + + if(scheck==1) warn("There were %d total q-records out where brute_near disagreed with cycle_for_near or find_near.",ncheck); + + return 0; + +} /* end of main for sunearqcsv */ + +/* -------------------------------------------------------------------------- */ +/* Specify compare function for qsort. */ + +int compSort1 (const void * q1, const void * q2) { + + struct QInfo* p1 = (struct QInfo*) q1; + struct QInfo* p2 = (struct QInfo*) q2; + + if(p1->dlots[locp] < p2->dlots[locp]) return (-1); + if(p1->dlots[locp] > p2->dlots[locp]) return (1); + + return (0); + +} +/* --------------------------------------------------------------------------------------------------- */ + +void connect_nodes (node *tree_nodes, unsigned long long tree_numc, + double **tree_dl, int tree_numd, int ihop) { + +/* This function connects the already allocated tree nodes. */ +/* */ +/* Input arguments: */ +/* tree_nodes The already fully allocated tree nodes. */ +/* tree_numc Number of nodes in tree_nodes. */ +/* tree_dl Pointers to first element of values of each dimension. That is, each coordinate is in */ +/* a separate contiguous array. These are the pointers to the start of those arrays. */ +/* Those arrays have size tree_numc. */ +/* tree_numd Number of pointers in tree_dl (i.e. number of dimensions). */ +/* ihop = 1 means process the points in dispersed order (approximately random). */ +/* The worst search times for kd trees occur when there are just a few long branches. */ +/* This happens in simple kd tree code if the points are ordered in increasing or */ +/* decreasing coordinates. Unfortunately, that is often the case for survey files. */ +/* The literature about kd trees explains how to create trees with branches of the same */ +/* depth (perfectly balanced trees). That is not done here. Instead, this option hops */ +/* through the points and loads them in more-or-less random order. This creates a */ +/* reasonably-balanced tree (unless you are very unlucky). */ +/* Note that seriously unbalanced trees still work, but they may be very slow. */ +/* = 0 process the points in the order they exist within their arrays. */ + + unsigned long long nrat = 0; + unsigned long long nd = 0; + unsigned long long n = 0; + unsigned long long m = 0; + + for(m=0; m0; nrat=nrat*0.6) { + if(nrat>6) nd = sqrt((double) nrat); + else nd = 0; + for (n=nd+nrat/2; nelem, which means the points */ +/* are actually loaded into the tree in a different order than the order in their arrays. */ + + unsigned long long m = 1; + int naxe = 0; + + for(m=1; melem] < tree_dl[naxe][now_node->elem]) + next_node = now_node->l; + else next_node = now_node->r; + } + + if(tree_dl[naxe][m_node->elem] < tree_dl[naxe][now_node->elem]) + now_node->l = m_node; + else now_node->r = m_node; + + } + + return; +} + +/* --------------------------------------------------------------------------------------------------- */ + +void find_in (node *tree_nodes, double **tree_dl, int tree_numd, + double *extent_min, double *extent_max, + unsigned long long *out_elem, unsigned long long *num_out) { + +/* This function finds all points between specified extents of the dimensions. */ +/* */ +/* Input arguments: */ +/* tree_nodes The already fully allocated and connected tree nodes. */ +/* tree_dl Pointers to first element of values of each dimension. That is, each coordinate is in */ +/* a separate contiguous array. These are the pointers to the start of those arrays. */ +/* Those arrays have size tree_numc. */ +/* tree_numd Number of pointers in tree_dl (i.e. number of dimensions). */ +/* extent_min Minimum value of this dimension to find points. Size tree_numd. */ +/* Greater OR EQUAL to this value is in range. */ +/* extent_max Maximum value of this dimension to find nearest point. Size tree_numd. */ +/* Strictly LESS than this value is in range. */ +/* */ +/* Output arguments: */ +/* out_elem The element numbers of the points in the tree_dl arrays. With m meaning the */ +/* coordinates of the point are tree_dl[n][m] where n is the dimension number. */ +/* num_out Is the number of points found within the extent ranges. */ +/* */ +/* Return:false means some kind of input argument error. */ +/* NOTE: inputting an impossible extent_min,extent_max range is not an error, you */ +/* just get 0 for num_out. */ + + int naxe = 1; + node *now_node; + + *num_out = 0; + now_node = tree_nodes; + find_in_rcur(tree_dl, tree_numd, now_node, naxe, extent_min, extent_max, out_elem, num_out); + + return; +} + +/* --------------------------------------------------------------------------------------------------- */ + +void find_in_rcur (double **tree_dl, int tree_numd, + node * now_node, int naxe, + double *extent_min, double *extent_max, + unsigned long long *out_elem, unsigned long long *num_out) { + +/* This function finds all points between specified extents of the dimensions. */ +/* */ +/* Input arguments: */ +/* tree_dl Pointers to first element of values of each dimension. That is, each coordinate is in */ +/* a separate contiguous array. These are the pointers to the start of those arrays. */ +/* Those arrays have size tree_numc. */ +/* tree_numd Number of pointers in tree_dl (i.e. number of dimensions). */ +/* extent_min Minimum value of this dimension to find points. Size tree_numd. */ +/* Greater OR EQUAL to this value is in range. */ +/* extent_max Maximum value of this dimension to find nearest point. Size tree_numd. */ +/* Strictly LESS than this value is in range. */ +/* */ +/* Input/Output arguments: */ +/* now_node The current node to start from. */ +/* naxe The current splitting dimension. */ +/* out_elem Accumulating element numbers of the points in the tree_dl arrays. With m meaning the */ +/* coordinates of the point are tree_dl[n][m] where n is the dimension number. */ +/* num_out Accumulating number of points found within the extent ranges. */ + + if(now_node==0) return; + + bool in = true; + int i = 0; + + for(i=0; ielem] < extent_min[i] || tree_dl[i][now_node->elem] >= extent_max[i]) { + in = false; + break; + } + } + if(in) { + out_elem[*num_out] = now_node->elem; + *num_out = *num_out + 1; + } + +/* Find more points. */ + + if(naxe==tree_numd) naxe = 0; + + if(tree_dl[naxe][now_node->elem] >= extent_min[naxe] && now_node->l) + find_in_rcur(tree_dl, tree_numd, now_node->l, naxe+1, extent_min, extent_max, out_elem, num_out); + + if(tree_dl[naxe][now_node->elem] < extent_max[naxe] && now_node->r) + find_in_rcur(tree_dl, tree_numd, now_node->r, naxe+1, extent_min, extent_max, out_elem, num_out); + + return; +} + +/* --------------------------------------------------------------------------------------------------- */ + +void find_near (node *tree_nodes, double **tree_dl, int *tree_nt, int tree_numd, + double *extent_min, double *extent_max, double *target, + unsigned long long *near_elem, double *near_dist, unsigned long long *num_found) { + +/* This function finds a nearest point between specified extents of the dimensions. */ +/* This function is usually slower for wide extents compared to function cycle_for_near. */ +/* */ +/* Input arguments: */ +/* tree_nodes The already fully allocated and connected tree nodes. */ +/* tree_dl Pointers to first element of values of each dimension. That is, each coordinate is in */ +/* a separate contiguous array. These are the pointers to the start of those arrays. */ +/* Those arrays have size tree_numc. */ +/* tree_nt Near type flags. 1 means standard pythagorean nearest. See note for what 0 means. */ +/* tree_numd Number of pointers in tree_dl and flags in tree_nt (i.e. number of dimensions). */ +/* extent_min Minimum value of this dimension to find nearest point. Size tree_numd. */ +/* Greater OR EQUAL to this value is in range. */ +/* extent_max Maximum value of this dimension to find nearest point. Size tree_numd. */ +/* Strictly LESS than this value is in range. */ +/* target Find the point nearest here considering extents and tree_nt. Size tree_numd. */ +/* */ +/* Output arguments: */ +/* near_elem The element number of a nearest point in the tree_dl arrays. For instance, a value of */ +/* m means the coordinates of point are tree_dl[n][m] where n is the dimension number. */ +/* near_dist The SQUARED distance between the target and nearest point (if num_found>0). */ +/* num_found 0 if no point is found within the extent ranges. */ +/* >0 is as many equally-near points as are found within the extent ranges. */ +/* The returned near_elem is the highest-numbered elem of all the nearest points. */ +/* */ +/* NOTE: A tree_nt value of 0 means the difference between the point coordinate and the target */ +/* coordinate are not added to the Pythagorean sum for the specified dimension(s). */ +/* So, the point distances are determined as if that dimension was not specified at all. */ +/* However, the extent_min,extent_max values for that dimension are still used. */ +/* The result is therefor the nearest point considering only the other dimensions, but */ +/* still restricted by the extent range of that dimension(s). */ + + int naxe = 1; + node *now_node; + + *num_found = 0; + *near_dist = DBL_MAX; + *near_elem = 0; + + now_node = tree_nodes; + find_near_rcur(tree_dl, tree_nt, tree_numd, now_node, naxe, + extent_min, extent_max, target, + near_elem, near_dist, num_found); + + return; +} + +/* --------------------------------------------------------------------------------------------------- */ + +void find_near_rcur (double **tree_dl, int *tree_nt, int tree_numd, node * now_node, int naxe, + double *extent_min, double *extent_max, double *target, + unsigned long long *near_elem, double *near_dist, unsigned long long *num_found) { + +/* This function finds a nearest point between specified extents of the dimensions. */ +/* This function is usually slower for wide extents compared to function cycle_for_near. */ +/* */ +/* Input arguments: */ +/* tree_dl Pointers to first element of values of each dimension. That is, each coordinate is in */ +/* a separate contiguous array. These are the pointers to the start of those arrays. */ +/* Those arrays have size tree_numc. */ +/* tree_nt Near type flags. 1 means standard pythagorean nearest. See note for what 0 means. */ +/* tree_numd Number of pointers in tree_dl and flags in tree_nt (i.e. number of dimensions). */ +/* extent_min Minimum value of this dimension to find nearest point. Size tree_numd. */ +/* Greater OR EQUAL to this value is in range. */ +/* extent_max Maximum value of this dimension to find nearest point. Size tree_numd. */ +/* Strictly LESS than this value is in range. */ +/* target Find the point nearest here considering extents and tree_nt. Size tree_numd. */ +/* */ +/* Input/Output arguments: */ +/* now_node The current node to start from. */ +/* naxe The current splitting dimension. */ +/* */ +/* Output arguments: */ +/* near_elem The element number of a nearest point in the tree_dl arrays. For instance, a value of */ +/* m means the coordinates of point are tree_dl[n][m] where n is the dimension number. */ +/* near_dist The SQUARED distance between the target and nearest point (if num_found>0). */ +/* num_found 0 if no point is found within the extent ranges. */ +/* >0 is as many equally-near points as are found within the extent ranges. */ +/* The returned near_elem is the highest-numbered elem of all the nearest points. */ +/* */ +/* NOTE: A tree_nt value of 0 means the difference between the point coordinate and the target */ +/* coordinate are not added to the Pythagorean sum for the specified dimension(s). */ +/* So, the point distances are determined as if that dimension was not specified at all. */ +/* However, the extent_min,extent_max values for that dimension are still used. */ +/* The result is therefor the nearest point considering only the other dimensions, but */ +/* still restricted by the extent range of that dimension(s). */ + + if(now_node==0) return; + + bool in = true; + int i = 0; + double rad = 0.; + + for(i=0; ielem] < extent_min[i] || tree_dl[i][now_node->elem] >= extent_max[i]) { + in = false; + break; + } + } + + if(in) { + rad = 0.; + for(i=0; ielem]) * (target[i]-tree_dl[i][now_node->elem]); + } + } + +/* The following set of ifs could be done differently. But I think coding it */ +/* this way saves cpu time by immediately rejecting the > cases using 1 if. */ +/* But who knows how modern optimizers will alter this. */ + + if(rad <= *near_dist) { + if(rad < *near_dist) { /* if distance is smaller, reset count to 1 */ + *num_found = 1; + *near_elem = now_node->elem; + *near_dist = rad; + } + else { /* so, distances are equal. Increment count, set output to higher elem.*/ + *num_found = *num_found + 1; + if(now_node->elem > *near_elem) *near_elem = now_node->elem; + } + } + + } + +/* Find more. */ + + if(naxe==tree_numd) naxe = 0; + + if(tree_dl[naxe][now_node->elem] >= extent_min[naxe] && now_node->l) + find_near_rcur(tree_dl, tree_nt, tree_numd, now_node->l, naxe+1, + extent_min, extent_max, target, + near_elem, near_dist, num_found); + + if(tree_dl[naxe][now_node->elem] < extent_max[naxe] && now_node->r) + find_near_rcur(tree_dl, tree_nt, tree_numd, now_node->r, naxe+1, + extent_min, extent_max, target, + near_elem, near_dist, num_found); + + return; +} + +/* --------------------------------------------------------------------------------------------------- */ + +void cycle_for_near(node *tree_nodes, double **tree_dl, int *tree_nt, int tree_numd, + double *extent_min, double *extent_max, double *target, + double init_rad, double rad_scal, + unsigned long long *near_elem, double *near_dist, + unsigned long long *num_found, int *ncycles) { + +/* This function finds a nearest point between specified extents of the dimensions. */ +/* This function is usually faster for wide extents compared to function find_near. */ +/* */ +/* Input arguments: */ +/* tree_nodes The already fully allocated and connected tree nodes. */ +/* tree_dl Pointers to first element of values of each dimension. That is, each coordinate is in */ +/* a separate contiguous array. These are the pointers to the start of those arrays. */ +/* Those arrays have size tree_numc. */ +/* tree_nt Near type flags. 1 means standard pythagorean nearest. See note for what 0 means. */ +/* tree_numd Number of pointers in tree_dl and flags in tree_nt (i.e. number of dimensions). */ +/* extent_min Minimum value of this dimension to find nearest point. Size tree_numd. */ +/* Greater OR EQUAL to this value is in range. */ +/* extent_max Maximum value of this dimension to find nearest point. Size tree_numd. */ +/* Strictly LESS than this value is in range. */ +/* target Find the point nearest here considering extents and tree_nt. Size tree_numd. */ +/* init_rad Initially search within radius init_rad. If no point within extents is found, */ +/* init_rad is multiplied by rad_scal and the seach is repeated. And so on. */ +/* rad_scal Scale value to apply to init_rad until a point is found within extents, or the */ +/* init_rad has been scaled to larger than the extents. */ +/* */ +/* Output arguments: */ +/* near_elem The element number of a nearest point in the tree_dl arrays. For instance, a value of */ +/* m means the coordinates of point are tree_dl[n][m] where n is the dimension number. */ +/* near_dist The SQUARED distance between the target and nearest point (if num_found>0). */ +/* num_found 0 if no point is found within the extent ranges. */ +/* >0 is as many equally-near points as are found within the extent ranges. */ +/* The returned near_elem is the highest-numbered elem of all the nearest points. */ +/* */ +/* ncycles Number of cycles to find nearest point. This is just for informational purposes. */ +/* This is the number of times that the rad_scal had to be applied before finding the */ +/* nearest point. This could help set init_rad and rad_scal for faster searches. */ +/* */ +/* NOTE: A tree_nt value of 0 means the difference between the point coordinate and the target */ +/* coordinate are not added to the Pythagorean sum for the specified dimension(s). */ +/* So, the point distances are determined as if that dimension was not specified at all. */ +/* However, the extent_min,extent_max values for that dimension are still used. */ +/* The result is therefor the nearest point considering only the other dimensions, but */ +/* still restricted by the extent range of that dimension(s). */ + + int naxe = 1; + node *now_node; + double now_rad = init_rad; + double loc_min[9]; + double loc_max[9]; + int nset = 0; + int i = 0; + int nbig = 0; + + for(i=0; i= extent_max[i]) { /* yes, >= is better than > here */ + loc_max[i] = extent_max[i]; + nbig++; + } + } + } + + *num_found = 0; + *near_dist = DBL_MAX; + *near_elem = 0; + + naxe = 1; + now_node = tree_nodes; + find_near_rcur(tree_dl, tree_nt, tree_numd, now_node, naxe, + loc_min, loc_max, target, + near_elem, near_dist, num_found); + + if(*num_found<1) { + if(nbig==2*tree_numd) return; /* None found. Are we outside the extents? */ + now_rad = now_rad * rad_scal; + goto CYCLE; + } + +/* Here, we need to consider the difference between a square and a circle. */ +/* The now_rad value is half the size of the square we just searched. So, if */ +/* current nearest point is in a circle with that radius, we are finished. */ +/* But, otherwise, there might be nearer points hiding in the area outside */ +/* the searched-square, but inside the circle. So increase the search size */ +/* to the CURRENT nearest point radius. On the next cycle, we will definitly */ +/* get the nearest point because the CURRENT point is definitly within the */ +/* square with the now_rad that we are now setting. So the CURRENT point here */ +/* will be among the points returned by the next cycle. And its radius will */ +/* definitly satisfy this condition because we explicitly made now_rad big */ +/* enough (but another point hiding in the square-circle area might sneak in).*/ + + if(*near_dist > now_rad*now_rad && nbig!=2*tree_numd) { + now_rad = sqrt(*near_dist) * 1.001; + goto CYCLE; + } + + return; +} + +/* --------------------------------------------------------------------------------------------------- */ + +void brute_near (double **tree_dl, unsigned long long tree_numc, int *tree_nt, int tree_numd, + double *extent_min, double *extent_max, double *target, + unsigned long long *near_elem, double *near_dist, unsigned long long *num_found) { + +/* This function finds a nearest point between specified extents of the dimensions. */ +/* This function is usually much slower than cycle_for_near and find_near. */ +/* The original purpose of this function was to confirm that cycle_for_near and find_near */ +/* worked as expected (modified kdtree searches are nothing to take for granted). */ +/* However, for a small number of points (tree_numc) this function will also be faster. */ +/* */ +/* Input arguments: */ +/* tree_dl Pointers to first element of values of each dimension. That is, each coordinate is in */ +/* a separate contiguous array. These are the pointers to the start of those arrays. */ +/* Those arrays have size tree_numc. */ +/* tree_numc Number of points. */ +/* tree_nt Near type flags. 1 means standard pythagorean nearest. See note for what 0 means. */ +/* tree_numd Number of pointers in tree_dl and flags in tree_nt (i.e. number of dimensions). */ +/* extent_min Minimum value of this dimension to find nearest point. Size tree_numd. */ +/* Greater OR EQUAL to this value is in range. */ +/* extent_max Maximum value of this dimension to find nearest point. Size tree_numd. */ +/* Strictly LESS than this value is in range. */ +/* target Find the point nearest here considering extents and tree_nt. Size tree_numd. */ +/* */ +/* Input/Output arguments: */ +/* */ +/* Output arguments: */ +/* near_elem The element number of a nearest point in the tree_dl arrays. For instance, a value of */ +/* m means the coordinates of point are tree_dl[n][m] where n is the dimension number. */ +/* near_dist The SQUARED distance between the target and nearest point (if num_found>0). */ +/* num_found 0 if no point is found within the extent ranges. */ +/* >0 is as many equally-near points as are found within the extent ranges. */ +/* The returned near_elem is the highest-numbered elem of all the nearest points. */ +/* */ +/* NOTE: A tree_nt value of 0 means the difference between the point coordinate and the target */ +/* coordinate are not added to the Pythagorean sum for the specified dimension(s). */ +/* So, the point distances are determined as if that dimension was not specified at all. */ +/* However, the extent_min,extent_max values for that dimension are still used. */ +/* The result is therefor the nearest point considering only the other dimensions, but */ +/* still restricted by the extent range of that dimension(s). */ + + *num_found = 0; + *near_dist = DBL_MAX; + *near_elem = 0; + + bool in = true; + int i = 0; + double rad = 0.; + + unsigned long long ibrute = 0; + + for(ibrute=0; ibrute= extent_max[i]) { + in = false; + break; + } + } + + if(in) { + rad = 0.; + for(i=0; i *near_elem) *near_elem = ibrute; + } + } + + } + + } /* end of for(ibrute=0; ibrute