*! version 1.3 jkc 31oct2014 /* bipolate provides bivariate interpolation and smooth surface fitting for values given at irregularly distributed points. Author: Joseph Canner Johns Hopkins University School of Medicine Department of Surgery Center for Surgical Trials and Outcomes Research jcanner1@jhmi.edu Version 1.0.0 May 19, 2014 Version 1.1 September 17, 2014 adds support for : (1) using a variable number of nearest neighbors for computing derivatives (old Akima method) (2) using the interpolation methods employed by twoway contour (3) varlist ordering consistent with twoway contour and ipolate Version 1.2 October 20, 2014 adds support for : (1) Thin-plate spline (Green & Silverman)--equivalent to Stata method (2) Thin-plate spline plus smoothing (3) Shepard--equivalent to Stata method (4) New option scheme ("interp" by default, "fill" implied by fillusing, convexhull available for all interp methods) Version 1.3 October 31, 2014 adds support for : (1) incorporating variability in z (i.e., multiple z values for a given (x,y) pair) */ program bipolate version 12 syntax varlist(min=3 max=3 numeric) [if] [in] [using] [, SAVing(string) METHod(string) CONTOURmethod(string) /// XGrid(numlist) XLevels(numlist integer >1 min=1 max=1) /// YGrid(numlist) YLevels(numlist >1 integer min=1 max=1) /// FILLusing(string) CONVEXhull COLLapse(string) /// KEEPMISSing NEAR(numlist >1 integer min=1 max=1) /// SMOOTH(real 0) REPS(integer 1) SEED(integer -1)] preserve tokenize `varlist' local z `1' local y `2' local x `3' if ("`using'"!="") { qui use `using', clear } if ("`if'"!="" | "`in'"!="") { qui keep `if' `in' } if ("`near'"!="") { if (`near'<2 | `near'>min(25,`=_N')) { di as red "Near must be be at least 2 and no larger than 25 or the number of observations, whichever is less" exit } } else { local near=-1 } if `reps'<1 { di as red "Reps must be at least 1" exit } // Implied option "mode" (interp or fill) local mode "interp" if ("`fillusing'"!="") { local mode "fill" } // Determine method if ("`method'"=="") { local method="akima" // The default } else { local method=lower("`method'") } // Check for contradictory options if ("`contourmethod'"!="" & "`method'"!="contour") { di as red "may not use option contourmethod() without method(contour)" exit } if (("`method'"=="contour" | "`contourmethod'"!="") & ("`xgrid'"!="" | "`ygrid'"!="" | "`xlevels'"!="" | "`ylevels'"!="")) { di as red "may not combine options method(contour)/contourmethod() and xgrid()/ygrid()/xlevels()/ylevels()" exit } if ("`fillusing'"!="" & ("`xgrid'"!="" | "`ygrid'"!="" | "`xlevels'"!="" | "`ylevels'"!="")) { di as red "may not combine options xgrid()/ygrid()/xlevels()/ylevels() and fillusing()" exit } if (("`xgrid'"!="" & "`xlevels'"!="") | ("`ygrid'"!="" & "`ylevels'"!="")) { di as red "may not combine option xgrid()/ygrid() and xlevels()/ylevels()" exit } if ("`fillusing'"!="" & "`convexhull'"!="") { di as red "Can't specify both fillusing and convexhull" exit } if `reps'>1 & "`collapse'"!="" { di as red "may not combine collapse() with reps()" exit } recast double `x' `y' `z' // If multiple z-values for a given (x,y) pair, choose a method to collapse and get rid of the duplicates (unless reps>1) if `reps'==1 { local collapse=lower("`collapse'") if ("`collapse'"=="") { local collapse="mean" } if ("`collapse'"!="none") { collapse (`collapse') `z'_`collapse'=`z', by(`x' `y') } else { rename `z' `z'_none } qui drop if `x'==. | `y'==. if ("`keepmissing'"=="") { qui drop if `z'_`collapse'==. } keep `x' `y' `z'_`collapse' } else { // do mean and sd // Need to add a check that ignores reps() if there are no (x,y) pairs with multiple z values collapse (mean) `z'_mean=`z' (sd) `z'_sd=`z', by(`x' `y') qui drop if `x'==. | `y'==. if ("`keepmissing'"=="") { qui drop if `z'_mean==. } keep `x' `y' `z'_mean `z'_sd } // Calculate xgrid from xlevels if ("`xlevels'"!="" & "`mode'"=="interp") { qui summ `x' local width=(`r(max)'-`r(min)')/(`xlevels'-1) numlist "`r(min)'(`width')`r(max)'" local xgrid "`r(numlist)'" //di "`xgrid'" } // Calculate ygrid from ylevels if ("`ylevels'"!="" & "`mode'"=="interp") { qui summ `y' local width=(`r(max)'-`r(min)')/(`ylevels'-1) numlist "`r(min)'(`width')`r(max)'" local ygrid "`r(numlist)'" //di "`ygrid'" } // Calculate xgrid from data if ("`xgrid'"=="" & "`xlevels'"=="" & "`mode'"=="interp") { qui levelsof `x', local(xgrid) } // Calculate ygrid from data if ("`ygrid'"=="" & "`ylevels'"=="" & "`mode'"=="interp") { qui levelsof `y', local(ygrid) } if ("`convexhull'"!="") { // Use convex hull method to generating a fillusing dataset // Already used one preserve command, use tempfile to preserve again tempfile preserve qui save "`preserve'", replace bip_cvxhull `y' `x', noreport nograph qui gen x1=`x' qui gen y1=min(_cvxh1l,_cvxh1r) // Keep just the hull vertices qui keep if !mi(x1) & !mi(y1) // Find approximate center of hull qui summ x1 local xmin=`r(min)' local xmax=`r(max)' local cx1=(`xmax'+`xmin')/2 qui summ y1 local ymin=`r(min)' local ymax=`r(max)' local cy1=(`ymax'+`ymin')/2 // Compute angles between center and each point and sort in angle order qui gen angle=atan2(y1-`cy1',x1-`cx1') sort angle // Find each points' next neighbor qui gen x2=x1[_n+1] qui gen y2=y1[_n+1] // Join the last point to the first point qui replace x2=x1[1] in l // l=last qui replace y2=y1[1] in l keep x1 x2 y1 y2 angle // For each potential point in the grid, determine if it is to the "left" of all sides of the hull. If so, skip it; if not, add it // See http://bbs.dartmouth.edu/~fangq/MATH/download/source/Determining%20if%20a%20point%20lies%20on%20the%20interior%20of%20a%20polygon.htm qui gen `x'=. qui gen `y'=. qui gen left=. qui count local nsides=`r(N)' // Loop through values of x and y implied by xgrid and ygrid foreach xx of numlist `xgrid' { foreach yy of numlist`ygrid' { //di "`xx' `yy'" // Determine if this potential point is to the "left" of all sides of the hull qui replace left=(`yy'-y1)*(x2-x1)-(`xx'-x1)*(y2-y1) // Count the number of points to the "left" qui count if left>=0 & !mi(left) // If all points are to the left, it is inside the hull, so add it to the list if `r(N)'==`nsides' { qui set obs `=_N+1' qui replace `x'=`xx' in `=_N' qui replace `y'=`yy' in `=_N' } } } // Drop the originall hull points since they are included in the above algorithm qui drop if !mi(angle) keep `x' `y' //scatter `y' `x' // Create a temporary fillusing file to use later tempfile fillusing qui save "`fillusing'", replace local mode "fill" use "`preserve'", clear } if ("`mode`'"=="fill") { qui merge 1:1 `x' `y' using `fillusing' drop _merge } if `reps'==1 { if ("`method'"=="contour") { // Stata method from twoway contour capture noisily call_contour `z'_`collapse' `y' `x', `contourmethod' if (_rc!=0) exit } else { capture noisily call_bivar `z'_`collapse' `y' `x', xgrid(`xgrid') ygrid(`ygrid') mode(`mode') `method' near(`near') smooth(`smooth') if (_rc!=0) exit } } else { // Preserve analysis file so that results can be accumulated tempfile preserve1 qui gen `z'_sim=. qui save `preserve1', replace tempfile results di "Reps:" _continue if `seed'>=0 { set seed `seed' } forval rep=1/`reps' { if mod(`rep',5)==0 { di "`rep'" _continue } else { di "." _continue } qui use `preserve1', clear qui replace `z'_sim=rnormal(`z'_mean,`z'_sd) qui replace `z'_sim=`z'_mean if inlist(`z'_sd,0,.) drop `z'_mean `z'_sd if ("`method'"=="contour") { // Stata method from twoway contour capture noisily call_contour `z'_sim `y' `x', `contourmethod' if (_rc!=0) exit } else { capture noisily call_bivar `z'_sim `y' `x', xgrid(`xgrid') ygrid(`ygrid') mode(`mode') `method' near(`near') smooth(`smooth') if (_rc!=0) exit } //desc //summ `z'_sim, det capture append using `results' qui save `results', replace } di _newline collapse (mean) `z'_mean=`z'_sim (sd) `z'_sd=`z'_sim, by(`x' `y') } if ("`saving'"!="") { save `saving' } end program call_bivar version 12 syntax varlist(min=3 max=3 numeric) , MODE(string) NEAR(integer) SMOOTH(real) /// [XGrid(numlist) YGrid(numlist) THINplatespline SHEPARD AKIMA] tokenize `varlist' local z `1' local y `2' local x `3' if ("`akima'"!="") { //di "Akima: `x' `y' `z' `xgrid' `ygrid' `mode' `near'" mata: bivar("`x'","`y'","`z'","`xgrid'","`ygrid'","`mode'",`near') } else if ("`thinplatespline'"!="") { //di "TPS: `x' `y' `z' `xgrid' `ygrid' `mode' `smooth'" mata: tps("`x'","`y'","`z'","`xgrid'","`ygrid'","`mode'",`smooth') } else if ("`shepard'"!="") { //di "Shep: `x' `y' `z' `xgrid' `ygrid' `mode' " mata: shep("`x'","`y'","`z'","`xgrid'","`ygrid'","`mode'") } end program call_contour version 12 syntax varlist(min=3 max=3 numeric) [, THINplatespline SHEPARD] tokenize `varlist' local z `1' local y `2' local x `3' if ("`thinplatespline'"!="" & "`shepard'"!="") { di as red "specify only one interpolation method for option contourmethod()" error 1 } if ("`shepard'"=="") { local contourmethod "tps" } else { local contourmethod "shepard" } mata: data=contour_interp("`z'","`y'","`x'","`contourmethod'") clear getmata (`1' `2' `3')=data, replace force end version 12 mata: mata set matastrict on real matrix contour_interp(string zvar, string yvar, string xvar, string method) { real matrix Z, data real scalar re, rn, dim, i, j Z = _contour_setup(zvar,yvar,xvar,newx=.,newy=. ,method) re = rows(newx) rn = rows(newy) dim = re*rn data = J(dim, 3, .) for(i=1; i<=re; i++) { for(j=1; j<=rn; j++) { pos = (i-1)*rn+j data[pos, .] = (Z[j,i] ,newy[j], newx[i]) } } return(data) } void bivar(string xvar, string yvar, string zvar, string xgrid, string ygrid, string mode, real scalar ncp) { real vector iwk, wk real vector x, y, z, xi, yi, zi, xd, yd, zd real matrix zij real scalar npt, ni, i, j, n x=st_data(.,xvar)' y=st_data(.,yvar)' z=st_data(.,zvar)' npt=st_nobs() if (mode=="interp") { real scalar nxi, nyi, nzi xi=strtoreal(tokens(xgrid)) nxi=cols(xi) yi=strtoreal(tokens(ygrid)) nyi=cols(yi) nzi=nxi zij=J(nzi,nyi,.) iwk=J(1,max((31,27+ncp))*npt+nxi*nyi,.) wk=J(1,6*npt,.) idsfft (1, ncp, npt, x, y, z, nxi, nyi, nzi, xi, yi, zij, iwk, wk) // Drop current observations and replace with newly interpolated observations st_dropobsin(.) st_addobs(nxi*nyi) n=0 for (i=1; i<=nxi; i++) { for (j=1; j<=nyi; j++) { n++ st_store(n,xvar, xi[i]) st_store(n,yvar, yi[j]) st_store(n,zvar, zij[i,j]) } } } else if (mode=="fill") { // Complete data set (original+filled) goes into xi,yi,zi ni=npt xi=x yi=y zi=J(1,ni,.) // Save points with non-missing z-values to be used for generating interpolation function norig=0 for (i=1; i<=npt; i++) { if (z[i]!=.) { norig++ xd=(xd,x[i]) yd=(yd,y[i]) zd=(zd,z[i]) } } npt=norig iwk=J(1,max((31,27+ncp))*npt+ni,.) wk=J(1,8*npt,.) idbvip ( 1, ncp, npt, xd, yd, zd, ni, xi, yi, zi, iwk, wk ) // Drop current observations and replace with newly filled-in observations (which will include the original observations) st_dropobsin(.) st_addobs(ni) st_store(.,xvar, xi') st_store(.,yvar, yi') st_store(.,zvar, zi') } return } void tps(string xvar, string yvar, string zvar, string xgrid, string ygrid, string mode, scalar alpha) { real vector iwk, wk real vector x, y, z, xi, yi, zi, xd, yd, zd real matrix zij real scalar n, ni, i, j, xx, yy x=st_data(.,xvar)' y=st_data(.,yvar)' z=st_data(.,zvar) n=st_nobs() if (mode=="interp") { real scalar nxi, nyi, nzi xi=strtoreal(tokens(xgrid)) nxi=cols(xi) yi=strtoreal(tokens(ygrid)) nyi=cols(yi) nzi=nxi zij=J(nzi,nyi,.) // TPS Code Ivec=J(1,n,1) t=x\y T=Ivec\x\y E=J(n,n,.) for (i=1; i<=n; i++) { for (j=1; j<=n; j++) { if (i==j) { E[i,j]=0 } else { E[i,j]=(1/(16*pi()))*(norm(t[.,i]-t[.,j])^2)*log(norm(t[.,i]-t[.,j])^2) } } } //"E=" //E mat0=J(3,3,0) vec0=J(3,1,0) A=(E+alpha*I(n,n),T'\T,mat0) B=(z\vec0) //"A= B=" // A //B X=lusolve(A,B) //"X=" //X d=X[1..n,.] a=X[n+1..n+3,.] //"d= a=" //d //a //"Tdelta" //T*d //"xi yi" //xi //yi for (xx=1; xx<=nxi; xx++) { for (yy=1; yy<=nyi; yy++) { /* "xx yy xy" xx yy */ xy=(xi[xx]\yi[yy]) //xy g1=0 for (i=1; i<=n; i++) { /* "i" i t[.,i] */ if (xy!=t[.,i]) { g1=g1+d[i]*(1/(16*pi()))*(norm(xy-t[.,i])^2)*log(norm(xy-t[.,i])^2) } } g2=a[1]+a[2]*xi[xx]+a[3]*yi[yy] zij[xx,yy]=g1+g2 } } // Drop current observations and replace with newly interpolated observations st_dropobsin(.) st_addobs(nxi*nyi) n=0 for (i=1; i<=nxi; i++) { for (j=1; j<=nyi; j++) { n++ st_store(n,xvar, xi[i]) st_store(n,yvar, yi[j]) st_store(n,zvar, zij[i,j]) } } } else if (mode=="fill") { // Complete data set (original+filled) goes into xi,yi,zi ni=n xi=x yi=y zi=J(1,ni,.) // Save points with non-missing z-values to be used for generating interpolation function norig=0 for (i=1; i<=n; i++) { if (z[i]!=.) { norig++ xd=(xd,x[i]) yd=(yd,y[i]) zd=(zd,z[i]) } } n=norig Ivec=J(1,n,1) t=xd\yd T=Ivec\xd\yd E=J(n,n,.) for (i=1; i<=n; i++) { for (j=1; j<=n; j++) { if (i==j) { E[i,j]=0 } else { E[i,j]=(1/(16*pi()))*(norm(t[.,i]-t[.,j])^2)*log(norm(t[.,i]-t[.,j])^2) } } } //"E=" //E mat0=J(3,3,0) vec0=J(3,1,0) A=(E+alpha*I(n,n),T'\T,mat0) //zd //vec0 B=(zd'\vec0) //"A= B=" //A //B X=lusolve(A,B) //"X=" //X d=X[1..n,.] a=X[n+1..n+3,.] //"d= a=" //d //a //"Tdelta" //T*d //"xi yi" //xi //yi for (i=1; i<=ni; i++) { xy=(xi[i]\yi[i]) g1=0 for (j=1; j<=n; j++) { if (xy!=t[.,j]) { g1=g1+d[j]*(1/(16*pi()))*(norm(xy-t[.,j])^2)*log(norm(xy-t[.,j])^2) } } g2=a[1]+a[2]*xi[i]+a[3]*yi[i] zi[i]=g1+g2 } // Drop current observations and replace with newly filled-in observations (which will include the original observations) st_dropobsin(.) st_addobs(ni) st_store(.,xvar, xi') st_store(.,yvar, yi') st_store(.,zvar, zi') } return } real scalar f(scalar y, scalar x) { return(1/(x^2+y^2)) } void shep(string xvar, string yvar, string zvar, string xgrid, string ygrid, string mode) { real vector iwk, wk real vector x, y, z, xi, yi, zi, xd, yd, zd real matrix zij real scalar npt, ni, i, j, n x=st_data(.,xvar)' y=st_data(.,yvar)' z=st_data(.,zvar) n=st_nobs() if (mode=="interp") { real scalar nxi, nyi, nzi xi=strtoreal(tokens(xgrid)) nxi=cols(xi) yi=strtoreal(tokens(ygrid)) nyi=cols(yi) nzi=nxi zij=J(nzi,nyi,.) // Shepard (from Stata documentation, G-2 page 200) for (xx=1; xx<=nxi; xx++) { for (yy=1; yy<=nyi; yy++) { denom=0 for (i=1; i<=n; i++) { denom=denom+f(yi[yy]-y[i],xi[xx]-x[i]) } numer=0 for (i=1; i<=n; i++) { numer=numer+z[i]*f(yi[yy]-y[i],xi[xx]-x[i]) } zij[xx,yy]=numer/denom } } // Drop current observations and replace with newly interpolated observations st_dropobsin(.) st_addobs(nxi*nyi) n=0 for (i=1; i<=nxi; i++) { for (j=1; j<=nyi; j++) { n++ st_store(n,xvar, xi[i]) st_store(n,yvar, yi[j]) st_store(n,zvar, zij[i,j]) } } } else if (mode=="fill") { // Complete data set (original+filled) goes into xi,yi,zi ni=n xi=x yi=y zi=J(1,ni,.) // Save points with non-missing z-values to be used for generating interpolation function norig=0 for (i=1; i<=n; i++) { if (z[i]!=.) { norig++ xd=(xd,x[i]) yd=(yd,y[i]) zd=(zd,z[i]) } } npt=norig for (i=1; i<=ni; i++) { denom=0 for (j=1; j<=n; j++) { denom=denom+f(yi[i]-yd[j],xi[i]-xd[j]) } numer=0 for (j=1; j<=n; j++) { numer=numer+zd[j]*f(yi[i]-yd[j],xi[i]-xd[j]) } zi[i]=numer/denom } // Drop current observations and replace with newly filled-in observations (which will include the original observations) st_dropobsin(.) st_addobs(ni) st_store(.,xvar, xi') st_store(.,yvar, yi') st_store(.,zvar, zi') } return } void idbvip ( real scalar md, real scalar ncp, real scalar ndp, real vector xd,real vector yd,real vector zd,real scalar nip,real vector xi,real vector yi,real vector zi, real vector iwk,real vector wk ) { /* ! !******************************************************************************* ! !! IDBVIP performs bivariate interpolation of irregular X, Y data. ! ! ! Discussion: ! ! The data points must be distinct and their projections in the ! X-Y plane must not be collinear, otherwise an error return ! occurs. ! ! Inadequate work space IWK and WK may cause incorrect results. ! ! Latest revision: ! ! January, 1985 ! ! Purpose: ! ! To provide bivariate interpolation and smooth surface fitting for ! values given at irregularly distributed points. ! ! The resulting interpolating function and its first-order partial ! derivatives are continuous. ! ! The method employed is local, i.e. a change in the data in one area ! of the plane does not affect the interpolating function except in ! that local area. This is advantageous over global interpolation ! methods. ! ! Also, the method gives exact results when all points lie in a plane. ! This is advantageous over other methods such as two-dimensional ! Fourier series interpolation. ! ! Usage: ! ! This package contains two user entries, IDBVIP and IDSFFT, both ! requiring input data to be given at points ! ( X(I), Y(I) ), I = 1,...,N. ! ! If the user desires the interpolated data to be output at grid ! points, i.e. at points ! ( XI(I), YI(J) ), I = 1,...,NXI, J=1,...,NYI, ! then IDSFFT should be used. This is useful for generating an ! interpolating surface. ! ! The other user entry point, IDBVIP, will produce interpolated ! values at scattered points ! ( XI(I), YI(I) ), i = 1,...,NIP. ! This is useful for filling in missing data points on a grid. ! ! History: ! ! The original version of BIVAR was written by Hiroshi Akima in ! August 1975 and rewritten by him in late 1976. It was incorporated ! into NCAR's public software libraries in January 1977. In August ! 1984 a new version of BIVAR, incorporating changes described in the ! Rocky Mountain Journal of Mathematics article cited below, was ! obtained from Dr Akima by Michael Pernice of NCAR's Scientific ! Computing Division, who evaluated it and made it available in February, ! 1985. ! ! Accuracy: ! ! Accurate to machine precision on the input data points. Accuracy at ! other points greatly depends on the input data. ! ! References: ! ! Hiroshi Akima, ! A Method of Bivariate Interpolation and Smooth Surface Fitting ! for Values Given at Irregularly Distributed Points, ! ACM Transactions on Mathematical Software, ! Volume 4, Number 2, June 1978. ! ! Hiroshi Akima, ! On Estimating Partial Derivatives for Bivariate Interpolation ! of Scattered Data, ! Rocky Mountain Journal of Mathematics, ! Volume 14, Number 1, Winter 1984. ! ! Method: ! ! The XY plane is divided into triangular cells, each cell having ! projections of three data points in the plane as its vertices, and ! a bivariate quintic polynomial in X and Y is fitted to each ! triangular cell. ! ! The coefficients in the fitted quintic polynomials are determined ! by continuity requirements and by estimates of partial derivatives ! at the vertices and along the edges of the triangles. The method ! described in the rocky mountain journal reference guarantees that ! the generated surface depends continuously on the triangulation. ! ! The resulting interpolating function is invariant under the following ! types of linear coordinate transformations: ! 1) a rotation of the XY coordinate system ! 2) linear scale transformation of the Z axis ! 3) tilting of the XY plane, i.e. new coordinates (u,v,w) given by ! u = x ! v = y ! w = z + a*x + b*y ! where a, b are arbitrary constants. ! ! complete details of the method are given in the reference publications. ! ! Parameters: ! ! Input, integer MD, mode of computation. MD must be 1, ! 2, or 3, else an error return occurs. ! ! 1: if this is the first call to this subroutine, or if the ! value of NDP has been changed from the previous call, or ! if the contents of the XD or YD arrays have been changed ! from the previous call. ! ! 2: if the values of NDP and the XD and YD arrays are unchanged ! from the previous call, but new values for XI, YI are being ! used. If MD = 2 and NDP has been changed since the previous ! call to IDBVIP, an error return occurs. ! ! 3: if the values of NDP, NIP, XD, YD, XI, YI are unchanged from ! the previous call, i.e. if the only change on input to idbvip is ! in the ZD array. If MD = 3 and NDP or NIP has been changed since ! the previous call to IDBVIP, an error return occurs. ! ! Between the call with MD = 2 or MD = 3 and the preceding call, the ! IWK and WK work arrays should not be disturbed. ! ! Input, integer NDP, the number of data points (must be 4 or ! greater, else an error return occurs). ! ! Input, real XD(NDP), YD(NDP), the X and Y coordinates of the data points. ! ! Input, real ZD(NDP), the data values at the data points. ! ! Input, integer NIP, the number of output points at which ! interpolation is to be performed (must be 1 or greater, else an ! error return occurs). ! ! Input, real XI(NIP), YI(NIP), the coordinates of the points at which ! interpolation is to be performed. ! ! Output, real ZI(NIP), the interpolated data values. ! ! Workspace, integer IWK(31*NDP+NIP). ! ! Workspace, real WK(8*NDP). ! */ // implicit none //real scalar ndp //real scalar nip //real scalar ap //real scalar bp //real scalar cp //real scalar dp real scalar iip //real scalar itipv //real scalar itpv //real vector iwk real scalar jwipl real scalar jwipt real scalar jwit real scalar jwit0 real scalar jwiwk real scalar jwiwl real scalar jwiwp real scalar jwwpd //real scalar md real scalar nl real scalar nt /* real scalar ntsc real scalar p00 real scalar p01 real scalar p02 real scalar p03 real scalar p04 real scalar p05 real scalar p10 real scalar p11 real scalar p12 real scalar p13 real scalar p14 real scalar p20 real scalar p21 real scalar p22 real scalar p23 real scalar p30 real scalar p31 real scalar p32 real scalar p40 real scalar p41 real scalar p50 */ //real vector wk //real scalar x0 //real vector xd //real vector xi //real scalar xs1 //real scalar xs2 //real scalar y0 //real vector yd //real vector yi //real scalar ys1 //real scalar ys2 //real vector zd //real vector zi real vector ipt, ipl, iwl, iwp, ipc real vector ngp, igp real vector wpd real scalar tzi real scalar it real vector tiwk // save /idlc/ // save /idpt/ // common /idlc/ itipv,xs1,xs2,ys1,ys2,ntsc(9) // common /idpt/ itpv,x0,y0,ap,bp,cp,dp, & // p00,p10,p20,p30,p40,p50,p01,p11,p21,p31,p41, & // p02,p12,p22,p32,p03,p13,p23,p04,p14,p05 external real scalar itipv,xs1,xs2,ys1,ys2 external real vector ntsc external real scalar itpv,x0,y0,ap,bp,cp,dp, p00,p10,p20,p30,p40,p50,p01,p11,p21,p31,p41, p02,p12,p22,p32,p03,p13,p23,p04,p14,p05 // Error check. if ( md < 1 | md > 3 ) { _error(1,"IDBVIP - Fatal error: Input parameter MD out of range.") return } if (ncp!=-1 & (ncp < 2 | ncp >= min((25,ndp))) ) { _error(6,"IDBVIP - Fatal error: Input parameter NCP out of range.") return } if ( ndp < 4 ) { _error(2,"IDBVIP - Fatal error: Input parameter NDP out of range.") return } if ( nip < 1 ) { _error(3,"IDBVIP - Fatal error: Input parameter NIP out of range.") return } if ( md == 1 ) { iwk[1] = ncp iwk[2] = ndp } else { if ( ncp != iwk[1] ) { _error(4,"IDBVIP - Fatal error: MD = 2 or 3 but NCP was changed since last call.") return } if ( ndp != iwk[2] ) { _error(5,"IDBVIP - Fatal error: MD = 2 or 3 but NDP was changed since last call.") return } } if ( md <= 2 ) { iwk[3] = nip } else { if ( nip < iwk[3] ) { _error (6,"IDBVIP - Fatal error: MD = 3 but NIP was changed since last call.") return } } // Allocation of storage areas in the IWK array. jwipt = 16 jwiwl = 6*ndp+1 jwiwk = jwiwl jwipl = 24*ndp+1 jwiwp = 30*ndp+1 jwipc = 27*ndp+1 jwit0 = max((31,27+ncp))*ndp jwwpd = 5*ndp+1 // Triangulate the XY plane. if ( md == 1 ) { ipt=iwk[|jwipt\.|] ipl=iwk[|jwipl\.|] iwl=iwk[|jwiwl\.|] iwp=iwk[|jwiwp\.|] // idtang ( ndp, xd, yd, nt, iwk[|jwipt\.|], nl, iwk[|jwipl\.|], iwk[|jwiwl\.|], iwk[|jwiwp\.|], wk ) idtang ( ndp, xd, yd, nt, ipt, nl, ipl, iwl, iwp, wk ) iwk[|jwipt\.|]=ipt iwk[|jwipl\.|]=ipl iwk[5] = nt iwk[6] = nl if ( nt == 0 ) { return } } else { nt = iwk[5] nl = iwk[6] } /* C DETERMINES NCP POINTS CLOSEST TO EACH DATA POINT. (FOR MD=1) 50 IF(MD0.GT.1) GO TO 60 CALL IDCLDP(NDP0,XD,YD,NCP0,IWK(JWIPC)) IF(IWK(JWIPC).EQ.0) RETURN */ if (md == 1 & ncp>=2) { idcldp(ndp, xd, yd, ncp, ipc) //iwk[|jwipc\.|]=ipc // For some reason this doesn't work here for (i=1; i<=cols(ipc); i++) { iwk[jwipc+i-1]=ipc[i] } if (iwk[jwipc] == 0) { return } } // Locate all points at which interpolation is to be performed. if ( md <= 2 ) { itipv = 0 jwit = jwit0 for (iip=1; iip<=nip; iip++) { jwit = jwit+1 ipt=iwk[|jwipt\.|] ipl=iwk[|jwipl\.|] tiwk=iwk[|jwiwk\.|] it=iwk[jwit] // idlctn ( ndp, xd, yd, nt, iwk[|jwipt\.|], nl, iwk[|jwipl\.|], xi[iip], yi[iip], iwk[jwit], iwk[|jwiwk\.|], wk ) idlctn ( ndp, xd, yd, nt, ipt, nl, ipl, xi[iip], yi[iip], it, tiwk, wk ) iwk[|jwiwk\.|]=tiwk iwk[jwit]=it } } // Estimate the partial derivatives at all data points. if (ncp==-1) { wpd=wk[|jwwpd\.|] // idpdrv ( ndp, xd, yd, zd, nt, iwk[|jwipt\.|], wk, wk[|jwwpd\.|] ) idpdrv ( ndp, xd, yd, zd, nt, ipt, wk, wpd ) wk[|jwwpd\.|]=wpd } else { ipc=iwk[|jwipc\.|] idpdrv_near ( ndp, xd, yd, zd, ncp, ipc, wk) iwk[|jwipc\.|]=ipc } // Interpolate the ZI values. itpv = 0 jwit = jwit0 for (iip=1; iip<=nip; iip++) { jwit = jwit + 1 //idptip ( ndp, xd, yd, zd, nt, iwk[jwipt], nl, iwk[jwipl], wk, iwk[jwit], xi[iip], yi[iip], zi[iip] ) tzi=. idptip ( ndp, xd, yd, zd, nt, ipt, nl, ipl, wk, iwk[jwit], xi[iip], yi[iip], tzi ) zi[iip]=tzi } return } real scalar spdt(real scalar u1, real scalar v1, real scalar u2, real scalar v2,real scalar u3,real scalar v3) { return((u1-u2)*(u3-u2)+(v1-v2)*(v3-v2)) } real scalar vpdt(real scalar u1, real scalar v1, real scalar u2, real scalar v2,real scalar u3,real scalar v3) { return((u1-u3)*(v2-v3)-(v1-v3)*(u2-u3)) } void idgrid ( real vector xd, real vector yd, real scalar nt, real vector ipt, real scalar nl, real vector ipl, real scalar nxi, real scalar nyi, real vector xi, real vector yi, real vector ngp, real vector igp ) { /* !******************************************************************************* ! !! IDGRID organizes grid points for surface fitting. ! ! ! Discussion: ! ! IDGRID sorts the points in ascending order of triangle numbers and ! of the border line segment number. ! ! Parameters: ! ! Input, real scalar XD(NDP), YD(NDP), the X and Y coordinates of the data ! points. ! ! Input, integer NT, the number of triangles. ! ! Input, integer IPT(3*NT), the indices of the triangle vertexes. ! ! Input, integer NL, the number of border line segments. ! ! Input, integer IPL(3*NL), containing the point numbers of the end points ! of the border line segments and their respective triangle numbers, ! ! Input, integer NXI, NYI, the number of grid points in the X and Y ! coordinates. ! ! Input, real scalar XI(NXI), YI(NYI), the coordinates of the grid points. ! ! Output, integer NGP(2*(NT+2*NL)) where the ! number of grid points that belong to each of the ! triangles or of the border line segments are to be stored. ! ! Output, integer IGP(NXI*NYI), where the grid point numbers are to be ! stored in ascending order of the triangle number and the border line ! segment number. ! */ // implicit none //real scalar nl //real scalar nt //real scalar nxi //real scalar nyi //real scalar igp(nxi*nyi) real scalar il0 real scalar il0t3 real scalar ilp1 real scalar ilp1t3 real scalar insd real scalar ip1 real scalar ip2 real scalar ip3 //real vector ipl(3*nl) //real vector ipt(3*nt) real scalar it0 real scalar it0t3 real scalar ixi real scalar iximn real scalar iximx real scalar iyi real scalar izi real scalar jigp0 real scalar jigp1 real scalar jigp1i real scalar jngp0 real scalar jngp1 real scalar l //real scalar ngp(2*(nt+2*nl)) real scalar ngp0 real scalar ngp1 real scalar nl0 real scalar nt0 real scalar nxinyi //real scalar spdt //real scalar u1 //real scalar u2 //real scalar u3 //real scalar v1 //real scalar v2 //real scalar v3 //real scalar vpdt real scalar x1 real scalar x2 real scalar x3 //real scalar xd(*) //real scalar xi(nxi) real scalar xii real scalar ximn real scalar ximx real scalar xmn real scalar xmx real scalar y1 real scalar y2 real scalar y3 //real scalar yd(*) //real scalar yi(nyi) real scalar yii real scalar yimn real scalar yimx real scalar ymn real scalar ymx /* ! ! Statement functions ! // Put these outside spdt(u1,v1,u2,v2,u3,v3) = (u1-u2)*(u3-u2)+(v1-v2)*(v3-v2) vpdt(u1,v1,u2,v2,u3,v3) = (u1-u3)*(v2-v3)-(v1-v3)*(u2-u3) */ // Preliminary processing nt0 = nt nl0 = nl nxinyi = nxi * nyi ximn = min (( xi[1], xi[nxi] )) ximx = max (( xi[1], xi[nxi] )) yimn = min (( yi[1], yi[nyi] )) yimx = max (( yi[1], yi[nyi] )) // Determine grid points inside the data area. jngp0 = 0 jngp1 = 2*(nt0+2*nl0)+1 jigp0 = 0 jigp1 = nxinyi + 1 for (it0=1; it0<=nt0; it0++) { ngp0 = 0 ngp1 = 0 it0t3 = it0*3 ip1 = ipt[it0t3-2] ip2 = ipt[it0t3-1] ip3 = ipt[it0t3] x1 = xd[ip1] y1 = yd[ip1] x2 = xd[ip2] y2 = yd[ip2] x3 = xd[ip3] y3 = yd[ip3] xmn = min (( x1, x2, x3 )) xmx = max (( x1, x2, x3 )) ymn = min (( y1, y2, y3 )) ymx = max (( y1, y2, y3 )) insd = 0 for (ixi=1; ixi<=nxi; ixi++) { if ( xi[ixi] < xmn | xi[ixi] > xmx ) { if ( insd == 0 ) { continue // cycle } iximx = ixi-1 goto L23 } if ( insd != 1 ) { insd = 1 iximn = ixi } } if ( insd == 0 ) { goto L38 } iximx = nxi L23: for (iyi = 1; iyi<=nyi; iyi++) { yii = yi[iyi] if ( yii < ymn | yii > ymx ) { goto L37 } for (ixi = iximn; ixi<=iximx; ixi++) { xii = xi[ixi] l = 0 if ( vpdt(x1,y1,x2,y2,xii,yii) < 0) goto L36 else if ( vpdt(x1,y1,x2,y2,xii,yii) == 0) goto L25 else if ( vpdt(x1,y1,x2,y2,xii,yii) > 0) goto L26 L25: l = 1 L26: if ( vpdt(x2,y2,x3,y3,xii,yii) < 0) goto L36 else if ( vpdt(x2,y2,x3,y3,xii,yii) == 0) goto L27 else if ( vpdt(x2,y2,x3,y3,xii,yii) > 0) goto L28 L27: l = 1 L28: if ( vpdt(x3,y3,x1,y1,xii,yii) < 0) goto L36 else if ( vpdt(x3,y3,x1,y1,xii,yii) == 0) goto L29 else if ( vpdt(x3,y3,x1,y1,xii,yii) > 0) goto L30 L29: l = 1 L30: izi = nxi*(iyi-1)+ixi if ( l == 1 ) goto L31 ngp0 = ngp0+1 jigp0 = jigp0+1 igp[jigp0] = izi goto L36 L31: for (jigp1i = jigp1; jigp1i<=nxinyi; jigp1i++) { if ( izi == igp[jigp1i] ) { goto L36 } } ngp1 = ngp1+1 jigp1 = jigp1-1 igp[jigp1] = izi L36: } L37: } L38: jngp0 = jngp0+1 ngp[jngp0] = ngp0 jngp1 = jngp1-1 ngp[jngp1] = ngp1 } // Determine grid points outside the data area. // in semi-infinite rectangular area. for (il0 = 1; il0<=nl0; il0++) { ngp0 = 0 ngp1 = 0 il0t3 = il0*3 ip1 = ipl[il0t3-2] ip2 = ipl[il0t3-1] x1 = xd[ip1] y1 = yd[ip1] x2 = xd[ip2] y2 = yd[ip2] xmn = ximn xmx = ximx ymn = yimn ymx = yimx if ( y2 >= y1 ) { xmn = min (( x1, x2 )) } if ( y2 <= y1 ) { xmx = max (( x1, x2 )) } if ( x2 <= x1 ) { ymn = min (( y1, y2 )) } if ( x2 >= x1 ) { ymx = max (( y1, y2 )) } insd = 0 for ( ixi = 1; ixi<=nxi; ixi++) { if ( xi[ixi] < xmn | xi[ixi] > xmx ) { if ( insd == 0 ) { goto L42 } iximx = ixi-1 goto L43 } if ( insd != 1 ) { insd = 1 iximn = ixi } L42: } if ( insd == 0 ) goto L58 iximx = nxi L43: for (iyi = 1; iyi<=nyi; iyi++) { yii = yi[iyi] if(yiiymx) goto L57 for (ixi = iximn; ixi<=iximx; ixi++) { xii = xi[ixi] l = 0 if ( vpdt(x1,y1,x2,y2,xii,yii) < 0) goto L46 else if ( vpdt(x1,y1,x2,y2,xii,yii) == 0) goto L45 else if ( vpdt(x1,y1,x2,y2,xii,yii) > 0) goto L56 L45: l = 1 L46: if ( spdt(x2,y2,x1,y1,xii,yii) < 0) goto L56 else if ( spdt(x2,y2,x1,y1,xii,yii) == 0) goto L47 else if ( spdt(x2,y2,x1,y1,xii,yii) > 0) goto L48 L47: l = 1 L48: if ( spdt(x1,y1,x2,y2,xii,yii) < 0) goto L56 else if ( spdt(x1,y1,x2,y2,xii,yii) == 0) goto L49 else if ( spdt(x1,y1,x2,y2,xii,yii) > 0) goto L50 L49: l = 1 L50: izi = nxi*(iyi-1)+ixi if(l==1) goto L51 ngp0 = ngp0+1 jigp0 = jigp0+1 igp[jigp0] = izi goto L56 L51: for (jigp1i = jigp1; jigp1i<=nxinyi; jigp1i++) { if(izi==igp[jigp1i]) goto L56 } L53: ngp1 = ngp1+1 jigp1 = jigp1-1 igp[jigp1] = izi L56: } L57: } L58: jngp0 = jngp0+1 ngp[jngp0] = ngp0 jngp1 = jngp1-1 ngp[jngp1] = ngp1 // In semi-infinite triangular area. L60: ngp0 = 0 ngp1 = 0 ilp1 = mod(il0,nl0)+1 ilp1t3 = ilp1*3 ip3 = ipl[ilp1t3-1] x3 = xd[ip3] y3 = yd[ip3] xmn = ximn xmx = ximx ymn = yimn ymx = yimx if(y3>=y2 & y2>=y1) xmn = x2 if(y3<=y2 & y2<=y1) xmx = x2 if(x3<=x2 & x2<=x1) ymn = y2 if(x3>=x2 & x2>=x1) ymx = y2 insd = 0 for (ixi = 1; ixi<=nxi; ixi++) { if ( xi[ixi] < xmn | xi[ixi] > xmx ) { if(insd==0) goto L62 iximx = ixi-1 goto L63 } if ( insd != 1 ) { insd = 1 iximn = ixi } L62: } if(insd==0) goto L78 iximx = nxi L63: for (iyi = 1; iyi<=nyi; iyi++) { yii = yi[iyi] if(yiiymx) goto L77 for (ixi = iximn; ixi<=iximx; ixi++) { xii = xi[ixi] l = 0 if ( spdt(x1,y1,x2,y2,xii,yii) < 0) goto L66 else if ( spdt(x1,y1,x2,y2,xii,yii) == 0) goto L65 else if ( spdt(x1,y1,x2,y2,xii,yii) > 0) goto L76 L65: l = 1 L66: if ( spdt(x3,y3,x2,y2,xii,yii) < 0) goto L70 else if ( spdt(x3,y3,x2,y2,xii,yii) == 0) goto L67 else if ( spdt(x3,y3,x2,y2,xii,yii) > 0) goto L76 L67: l = 1 L70: izi = nxi*(iyi-1)+ixi if ( l != 1 ) { ngp0 = ngp0+1 jigp0 = jigp0+1 igp[jigp0] = izi goto L76 } for (jigp1i = jigp1; jigp1i<=nxinyi; jigp1i++) { if(izi==igp[jigp1i]) goto L76 } ngp1 = ngp1+1 jigp1 = jigp1-1 igp[jigp1] = izi L76: } L77: } L78: jngp0 = jngp0+1 ngp[jngp0] = ngp0 jngp1 = jngp1-1 ngp[jngp1] = ngp1 } return } void idlctn (real scalar ndp, real vector xd, real vector yd, real scalar nt, real vector ipt, real scalar nl, real vector ipl, real scalar xii, real scalar yii, real scalar iti, real vector iwk, real vector wk ) { /* ! !******************************************************************************* ! !! IDLCTN finds the triangle that contains a point. ! ! ! Discusstion: ! ! IDLCTN determines what triangle a given point (XII, YII) belongs to. ! When the given point does not lie inside the data area, IDLCTN ! determines the border line segment when the point lies in an outside ! rectangular area, and two border line segments when the point ! lies in an outside triangular area. ! ! Parameters: ! ! Input, integer NDP, the number of data points. ! ! Input, real scalar XD(NDP), YD(NDP), the X and Y coordinates of the data. ! ! Input, integer NT, the number of triangles. ! ! Input, integer IPT(3*NT), the point numbers of the vertexes of ! the triangles, ! ! Input, integer NL, the number of border line segments. ! ! Input, integer IPL(3*NL), the point numbers of the end points of ! the border line segments and their respective triangle numbers. ! ! Input, real scalar XII, YII, the coordinates of the point to be located. ! ! Output, integer ITI, the triangle number, when the point is inside the ! data area, or two border line segment numbers, il1 and il2, ! coded to il1*(nt+nl)+il2, when the point is outside the data area. ! ! Workspace, integer IWK(18*NDP). ! ! Workspace, real scalar WK(8*NDP). ! */ // implicit none //real scalar ndp //real scalar nl //real scalar nt real scalar i1 real scalar i2 real scalar i3 real scalar idp real vector idsc real scalar il1 real scalar il1t3 real scalar il2 real scalar ip1 real scalar ip2 real scalar ip3 //real scalar ipl(3*nl) //real scalar ipt(3*nt) real scalar isc real scalar it0 real scalar it0t3 //real scalar iti //real scalar itipv real scalar itsc //real vector iwk real scalar jiwk real scalar jwk real scalar nl0 real scalar nt0 real scalar ntl //real vector ntsc real scalar ntsci //real scalar spdt //real scalar u1 //real scalar u2 //real scalar u3 //real scalar v1 //real scalar v2 //real scalar v3 //real scalar vpdt //real vector wk real scalar x0 real scalar x1 real scalar x2 real scalar x3 //real scalar xd(ndp) //real scalar xii real scalar xmn real scalar xmx //real scalar xs1 //real scalar xs2 real scalar y0 real scalar y1 real scalar y2 real scalar y3 //real scalar yd(ndp) //real scalar yii real scalar ymn real scalar ymx //real scalar ys1 //real scalar ys2 // save /idlc/ // common /idlc/ itipv,xs1,xs2,ys1,ys2,ntsc(9) external real scalar itipv,xs1,xs2,ys1,ys2 external real vector ntsc /* ! ! Statement functions ! spdt(u1,v1,u2,v2,u3,v3) = (u1-u2)*(u3-u2)+(v1-v2)*(v3-v2) vpdt(u1,v1,u2,v2,u3,v3) = (u1-u3)*(v2-v3)-(v1-v3)*(u2-u3) */ // Preliminary processing nt0 = nt nl0 = nl ntl = nt0+nl0 x0 = xii y0 = yii // Processing for a new set of data points if ( itipv!=0) goto L30 // Divide the x-y plane into nine rectangular sections. xmn = xd[1] xmx = xd[1] ymn = yd[1] ymx = yd[1] for (idp = 2; idp<=ndp; idp++) { xmn = min (( xd[idp], xmn )) xmx = max (( xd[idp], xmx )) ymn = min (( yd[idp], ymn )) ymx = max (( yd[idp], ymx )) } xs1 = ( xmn + xmn + xmx ) / 3.0E+00 xs2 = ( xmn + xmx + xmx ) / 3.0E+00 ys1 = ( ymn + ymn + ymx ) / 3.0E+00 ys2 = ( ymn + ymx + ymx ) / 3.0E+00 // Determine and store in the iwk array, triangle numbers of // the triangles associated with each of the nine sections. ntsc=J(1,9,0) idsc=J(1,9,0) it0t3 = 0 jwk = 0 for (it0 = 1; it0<=nt0; it0++) { it0t3 = it0t3+3 i1 = ipt[it0t3-2] i2 = ipt[it0t3-1] i3 = ipt[it0t3] xmn = min (( xd[i1], xd[i2], xd[i3] )) xmx = max (( xd[i1], xd[i2], xd[i3] )) ymn = min (( yd[i1], yd[i2], yd[i3] )) ymx = max (( yd[i1], yd[i2], yd[i3] )) if ( ymn <= ys1 ) { if(xmn<=xs1) idsc[1] = 1 if(xmx>=xs1 & xmn<=xs2) idsc[2] = 1 if(xmx>=xs2) idsc[3] = 1 } if ( ymx >= ys1 & ymn <= ys2 ) { if(xmn<=xs1) idsc[4] = 1 if(xmx>=xs1 & xmn<=xs2) idsc[5] = 1 if(xmx>=xs2) idsc[6] = 1 } if(ymx=xs1 & xmn<=xs2) idsc[8] = 1 if(xmx>=xs2) idsc[9] = 1 L25: for (isc = 1; isc<=9; isc++) { if ( idsc[isc] != 0 ) { jiwk = 9*ntsc[isc]+isc iwk[jiwk] = it0 ntsc[isc] = ntsc[isc]+1 idsc[isc] = 0 } } // Store in the wk array the minimum and maximum of the X and // Y coordinate values for each of the triangle. jwk = jwk+4 wk[jwk-3] = xmn wk[jwk-2] = xmx wk[jwk-1] = ymn wk[jwk] = ymx } goto L60 // Check if in the same triangle as previous. L30: it0 = itipv if(it0>nt0) goto L40 it0t3 = it0*3 ip1 = ipt[it0t3-2] x1 = xd[ip1] y1 = yd[ip1] ip2 = ipt[it0t3-1] x2 = xd[ip2] y2 = yd[ip2] if(vpdt(x1,y1,x2,y2,x0,y0) < 0.0E+00 ) goto L60 ip3 = ipt[it0t3] x3 = xd[ip3] y3 = yd[ip3] if(vpdt(x2,y2,x3,y3,x0,y0) < 0.0E+00 ) goto L60 if(vpdt(x3,y3,x1,y1,x0,y0) < 0.0E+00 ) goto L60 iti = it0 itipv = it0 return // Check if on the same border line segment. L40: il1 = trunc(it0 / ntl) il2 = it0-il1*ntl il1t3 = il1*3 ip1 = ipl[il1t3-2] x1 = xd[ip1] y1 = yd[ip1] ip2 = ipl[il1t3-1] x2 = xd[ip2] y2 = yd[ip2] if(il2!=il1) goto L50 if(spdt(x1,y1,x2,y2,x0,y0) < 0.0E+00 ) goto L60 if(spdt(x2,y2,x1,y1,x0,y0) < 0.0E+00 ) goto L60 if(vpdt(x1,y1,x2,y2,x0,y0) > 0.0E+00 ) goto L60 iti = it0 itipv = it0 return // Check if between the same two border line segments. L50: if(spdt(x1,y1,x2,y2,x0,y0) > 0.0E+00 ) goto L60 ip3 = ipl[3*il2-1] x3 = xd[ip3] y3 = yd[ip3] if ( spdt(x3,y3,x2,y2,x0,y0) <= 0.0E+00 ) { iti = it0 itipv = it0 return } // Locate inside the data area. // Determine the section in which the point in question lies. L60: isc = 1 if ( x0 >= xs1 ) { isc = isc+1 } if ( x0 >= xs2 ) { isc = isc+1 } if ( y0 >= ys1 ) { isc = isc+3 } if ( y0 >= ys2 ) { isc = isc+3 } // Search through the triangles associated with the section. ntsci = ntsc[isc] if(ntsci<=0) goto L70 jiwk = -9+isc for (itsc = 1; itsc<=ntsci; itsc++) { jiwk = jiwk+9 it0 = iwk[jiwk] jwk = it0*4 if(x0wk[jwk-2]) goto L61 if(y0wk[jwk]) goto L61 it0t3 = it0*3 ip1 = ipt[it0t3-2] x1 = xd[ip1] y1 = yd[ip1] ip2 = ipt[it0t3-1] x2 = xd[ip2] y2 = yd[ip2] if(vpdt(x1,y1,x2,y2,x0,y0)<0.0E+00 ) goto L61 ip3 = ipt[it0t3] x3 = xd[ip3] y3 = yd[ip3] if ( vpdt(x2,y2,x3,y3,x0,y0) >= 0.0E+00 ) { if ( vpdt(x3,y3,x1,y1,x0,y0) >= 0.0E+00 ) { iti = it0 itipv = it0 return } } L61: } // Locate outside the data area. L70: for (il1 = 1; il1<=nl0; il1++) { il1t3 = il1*3 ip1 = ipl[il1t3-2] x1 = xd[ip1] y1 = yd[ip1] ip2 = ipl[il1t3-1] x2 = xd[ip2] y2 = yd[ip2] if(spdt(x2,y2,x1,y1,x0,y0)<0.0E+00 ) goto L72 if(spdt(x1,y1,x2,y2,x0,y0)<0.0E+00 ) goto L71 if(vpdt(x1,y1,x2,y2,x0,y0)>0.0E+00 ) goto L72 il2 = il1 goto L75 L71: il2 = mod(il1,nl0)+1 ip3 = ipl[3*il2-1] x3 = xd[ip3] y3 = yd[ip3] if(spdt(x3,y3,x2,y2,x0,y0)<=0.0E+00 ) goto L75 L72: } it0 = 1 iti = it0 itipv = it0 return L75: it0 = il1*ntl+il2 iti = it0 itipv = it0 return } real scalar dsqf(real scalar u1, real scalar v1, real scalar u2, real scalar v2) { return((u2-u1)^2+(v2-v1)^2) } void idpdrv ( real scalar ndp, real vector xd, real vector yd, real vector zd, real scalar nt, real vector ipt, real vector pd, real vector wk ) { /* ! !******************************************************************************* ! !! IDPDRV estimates first and second partial derivatives at data points. ! ! ! Parameters: ! ! Input, integer NDP, the number of data points. ! ! Input, real scalar XD(NDP), YD(NDP), the X and Y coordinates of the data. ! ! Input, real scalar ZD(NDP), the data values. ! ! Input, integer NT, the number of triangles. ! ! Input, integer IPT(3*NT), the point numbers of the vertexes of the ! triangles. ! ! Output, real scalar PD(5*NDP), the estimated zx, zy, zxx, zxy, and zyy values ! at the ith data point are to be stored as the (5*i-4)th, (5*i-3)rd, ! (5*i-2)nd, (5*i-1)st and (5*i)th elements, respectively, where i = ! 1, 2, ..., ndp. ! ! Workspace, real scalar WK(NDP). ! */ // implicit none //integer ndp //integer nt real scalar d12 real scalar d23 real scalar d31 real scalar dx1 real scalar dx2 real scalar dy1 real scalar dy2 real scalar dz1 real scalar dz2 real scalar dzx1 real scalar dzx2 real scalar dzy1 real scalar dzy2 real scalar epsln epsln = 1.0E-06 real scalar idp //real scalar ipt(3*nt) real vector ipti ipti=J(1,3,.) real scalar it real scalar iv real scalar jpd real scalar jpd0 real scalar jpdmx real scalar jpt real scalar jpt0 real scalar nt0 //real scalar pd(5*ndp) real scalar vpx real scalar vpxx real scalar vpxy real scalar vpy real scalar vpyx real scalar vpyy real scalar vpz real scalar vpzmn real vector w1 w1=J(1,3,.) real scalar w2 w2=J(1,3,.) real scalar wi //real scalar wk(ndp) //real scalar xd(ndp) real scalar xv xv=J(1,3,.) //real scalar yd(ndp) real scalar yv yv=J(1,3,.) //real scalar zd(ndp) real scalar zv zv=J(1,3,.) real scalar zxv zxv=J(1,3,.) real scalar zyv zyv=J(1,3,.) // Preliminary processing. nt0 = nt // Clear the PD array. jpdmx = 5*ndp pd[1..jpdmx]=J(1,jpdmx,0.0E+00) wk[1..ndp]=J(1,ndp,0.0E+00) // Estimate ZX and ZY. for (it = 1; it<=nt0; it++) { jpt0 = 3*(it-1) for (iv = 1; iv<=3; iv++) { jpt = jpt0+iv idp = ipt[jpt] ipti[iv] = idp xv[iv] = xd[idp] yv[iv] = yd[idp] zv[iv] = zd[idp] } // Compute Vector Product of P1P2 and P1P3 dx1 = xv[2]-xv[1] dy1 = yv[2]-yv[1] dz1 = zv[2]-zv[1] dx2 = xv[3]-xv[1] dy2 = yv[3]-yv[1] dz2 = zv[3]-zv[1] vpx = dy1*dz2-dz1*dy2 vpy = dz1*dx2-dx1*dz2 vpz = dx1*dy2-dy1*dx2 vpzmn = abs(dx1*dx2+dy1*dy2)*epsln // Adjustment based on Rocky Mountain paper if ( abs(vpz) > vpzmn ) { // Compute distance between each pair of points in the triangle d12 = sqrt((xv[2]-xv[1])^2+(yv[2]-yv[1])^2) d23 = sqrt((xv[3]-xv[2])^2+(yv[3]-yv[2])^2) d31 = sqrt((xv[1]-xv[3])^2+(yv[1]-yv[3])^2) // Compute weights (reciprocal of the product of the paired distances) w1[1] = 1.0E+00 / (d31*d12) w1[2] = 1.0E+00 / (d12*d23) w1[3] = 1.0E+00 / (d23*d31) w2[1] = vpz*w1[1] w2[2] = vpz*w1[2] w2[3] = vpz*w1[3] // Sum up 3 vectors for (iv = 1; iv<=3; iv++) { idp = ipti[iv] jpd0 = 5*(idp-1) wi = (w1[iv]^2)*w2[iv] pd[jpd0+1] = pd[jpd0+1]+vpx*wi pd[jpd0+2] = pd[jpd0+2]+vpy*wi wk[idp] = wk[idp]+vpz*wi } } } for (idp = 1; idp<=ndp; idp++) { jpd0 = 5*(idp-1) pd[jpd0+1] = -pd[jpd0+1]/wk[idp] pd[jpd0+2] = -pd[jpd0+2]/wk[idp] } // Estimate ZXX, ZXY, and ZYY. for (it = 1; it<=nt0; it++) { jpt0 = 3*(it-1) for (iv = 1; iv<=3; iv++) { jpt = jpt0+iv idp = ipt[jpt] ipti[iv] = idp xv[iv] = xd[idp] yv[iv] = yd[idp] jpd0 = 5*(idp-1) zxv[iv] = pd[jpd0+1] zyv[iv] = pd[jpd0+2] } dx1 = xv[2]-xv[1] dy1 = yv[2]-yv[1] dzx1 = zxv[2]-zxv[1] dzy1 = zyv[2]-zyv[1] dx2 = xv[3]-xv[1] dy2 = yv[3]-yv[1] dzx2 = zxv[3]-zxv[1] dzy2 = zyv[3]-zyv[1] vpxx = dy1*dzx2-dzx1*dy2 vpxy = dzx1*dx2-dx1*dzx2 vpyx = dy1*dzy2-dzy1*dy2 vpyy = dzy1*dx2-dx1*dzy2 vpz = dx1*dy2-dy1*dx2 vpzmn = abs(dx1*dx2+dy1*dy2)*epsln if ( abs(vpz) > vpzmn ) { d12 = sqrt((xv[2]-xv[1])^2+(yv[2]-yv[1])^2) d23 = sqrt((xv[3]-xv[2])^2+(yv[3]-yv[2])^2) d31 = sqrt((xv[1]-xv[3])^2+(yv[1]-yv[3])^2) w1[1] = 1.0E+00 /(d31*d12) w1[2] = 1.0E+00 /(d12*d23) w1[3] = 1.0E+00 /(d23*d31) w2[1] = vpz*w1[1] w2[2] = vpz*w1[2] w2[3] = vpz*w1[3] for (iv = 1; iv<=3; iv++) { idp = ipti[iv] jpd0 = 5*(idp-1) wi = (w1[iv]^2)*w2[iv] pd[jpd0+3] = pd[jpd0+3]+vpxx*wi pd[jpd0+4] = pd[jpd0+4]+(vpxy+vpyx)*wi pd[jpd0+5] = pd[jpd0+5]+vpyy*wi } } } for (idp = 1; idp<=ndp; idp++) { jpd0 = 5*(idp-1) pd[jpd0+3] = -pd[jpd0+3]/wk[idp] pd[jpd0+4] = -pd[jpd0+4]/(2.0*wk[idp]) pd[jpd0+5] = -pd[jpd0+5]/wk[idp] } return } void idcldp(real scalar NDP, real vector XD, real vector YD, real scalar NCP, real vector IPC) { // Subroutine required by old version of IDPDRV /* SUBROUTINE IDCLDP(NDP,XD,YD,NCP,IPC) ID002720 C THIS SUBROUTINE SELECTS SEVERAL DATA POINTS THAT ARE CLOSEST C TO EACH OF THE DATA POINT. C THE INPUT PARAMETERS ARE C NDP = NUMBER OF DATA POINTS, C XD,YD = ARRAYS OF DIMENSION NDP CONTAINING THE X AND Y C COORDINATES OF THE DATA POINTS, C NCP = NUMBER OF DATA POINTS CLOSEST TO EACH DATA C POINTS. C THE OUTPUT PARAMETER IS C IPC = INTEGER ARRAY OF DIMENSION NCP*NDP, WHERE THE C POINT NUMBERS OF NCP DATA POINTS CLOSEST TO C EACH OF THE NDP DATA POINTS ARE TO BE STORED. C THIS SUBROUTINE ARBITRARILY SETS A RESTRICTION THAT NCP MUST C NOT EXCEED 25. C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS, C THEREFORE, SYSTEM DEPENDENT. C DECLARATION STATEMENTS DIMENSION XD(100),YD(100),IPC(400) DIMENSION DSQ0(25),IPC0(25) DATA NCPMX/25/, LUN/6/ C STATEMENT FUNCTION DSQF(U1,V1,U2,V2)=(U2-U1)**2+(V2-V1)**2 */ real scalar NCPMX NCPMX=25 real vector DSQ0, IPC0 DSQ0=J(1,NCPMX,.) IPC0=J(1,NCPMX,.) IPC=J(1,NCP*NDP,.) //C PRELIMINARY PROCESSING L10: NDP0=NDP NCP0=NCP if((NDP0<2) | (NCP0<1) | (NCP0>NCPMX) | (NCP0>=NDP0)) { _error(2090,"IDCLPD: IMPROPER INPUT PARAMETER VALUE(S)") return } //C CALCULATION L20: for (IP1=1; IP1<=NDP0; IP1++) { //C - SELECTS NCP POINTS. X1=XD[IP1] Y1=YD[IP1] J1=0 DSQMX=0.0 for (IP2=1; IP2<=NDP0; IP2++) { if(IP2 == IP1) goto L22 DSQI=dsqf(X1,Y1,XD[IP2],YD[IP2]) J1=J1+1 DSQ0[J1]=DSQI IPC0[J1]=IP2 if(DSQI <= DSQMX) goto L21 DSQMX=DSQI JMX=J1 L21: if(J1 >= NCP0) goto L23 L22: } L23: IP2MN=IP2+1 if(IP2MN > NDP0) goto L30 for (IP2=IP2MN; IP2<=NDP0; IP2++) { if(IP2 == IP1) goto L25 DSQI=dsqf(X1,Y1,XD[IP2],YD[IP2]) if(DSQI >= DSQMX) goto L25 DSQ0[JMX]=DSQI IPC0[JMX]=IP2 DSQMX=0.0 for (J1=1; J1<=NCP0; J1++) { if(DSQ0[J1] <= DSQMX) goto L24 DSQMX=DSQ0[J1] JMX=J1 L24: } L25: } //C - CHECKS IF ALL THE NCP+1 POINTS ARE COLLINEAR. L30: IP2=IPC0[1] DX12=XD[IP2]-X1 DY12=YD[IP2]-Y1 for (J3=2;J3<=NCP0; J3++) { IP3=IPC0[J3] DX13=XD[IP3]-X1 DY13=YD[IP3]-Y1 if((DY13*DX12-DX13*DY12) != 0.0) goto L50 } //C - SEARCHES FOR THE CLOSEST NONCOLLINEAR POINT. L40: NCLPT=0 for (IP3=1;IP3<=NDP0; IP3++) { if(IP3 == IP1) goto L43 for (J4=1;J4<=NCP0;J4++) { if(IP3 == IPC0[J4]) goto L43 } DX13=XD[IP3]-X1 DY13=YD[IP3]-Y1 if((DY13*DX12-DX13*DY12) == 0.0) goto L43 DSQI=dsqf(X1,Y1,XD[IP3],YD[IP3]) if(NCLPT == 0) goto L42 if(DSQI >= DSQMN) goto L43 L42: NCLPT=1 DSQMN=DSQI IP3MN=IP3 L43: } if(NCLPT == 0) { _error(2091,"IDCLPD: ALL COLLINEAR DATA POINTS") return } DSQMX=DSQMN IPC0[JMX]=IP3MN //C - REPLACES THE LOCAL ARRAY FOR THE OUTPUT ARRAY. L50: J1=(IP1-1)*NCP0 for (J2=1;J2<=NCP0;J2++) { J1=J1+1 IPC[J1]=IPC0[J2] } } /* C ERROR EXIT 90 WRITE (LUN,2090) GO TO 92 91 WRITE (LUN,2091) 92 WRITE (LUN,2092) NDP0,NCP0 IPC(1)=0 RETURN */ return } void idpdrv_near ( real scalar NDP, real vector XD, real vector YD, real vector ZD, real scalar NCP, real vector IPC, real vector PD ) { // Old version of IDPDRV that allows nearest neighbor derivatives /* C THIS SUBROUTINE ESTIMATES PARTIAL DERIVATIVES OF THE FIRST AND C SECOND ORDER AT THE DATA POINTS. C THE INPUT PARAMETERS ARE C NDP = NUMBER OF DATA POINTS, C XD,YD,ZD = ARRAYS OF DIMENSION NDP CONTAINING THE X, C Y, AND Z COORDINATES OF THE DATA POINTS, C NCP = NUMBER OF ADDITIONAL DATA POINTS USED FOR ESTI- C MATING PARTIAL DERIVATIVES AT EACH DATA POINT, C IPC = INTEGER ARRAY OF DIMENSION NCP*NDP CONTAINING C THE POINT NUMBERS OF NCP DATA POINTS CLOSEST TO C EACH OF THE NDP DATA POINTS. C THE OUTPUT PARAMETER IS C PD = ARRAY OF DIMENSION 5*NDP, WHERE THE ESTIMATED C ZX, ZY, ZXX, ZXY, AND ZYY VALUES AT THE DATA C POINTS ARE TO BE STORED. */ //C DECLARATION STATEMENTS // DIMENSION XD(100),YD(100),ZD(100),IPC(400),PD(500) real scalar NMX,NMY,NMZ,NMXX,NMXY,NMYX,NMYY PD=J(1,5*NDP,.) //C PRELIMINARY PROCESSING L10: NDP0=NDP NCP0=NCP NCPM1=NCP0-1 //C ESTIMATION OF ZX AND ZY L20: for (IP0=1; IP0<=NDP0; IP0++) { X0=XD[IP0] Y0=YD[IP0] Z0=ZD[IP0] NMX=0.0 NMY=0.0 NMZ=0.0 JIPC0=NCP0*(IP0-1) for (IC1=1; IC1<=NCPM1; IC1++) { JIPC=JIPC0+IC1 IPI=IPC[JIPC] DX1=XD[IPI]-X0 DY1=YD[IPI]-Y0 DZ1=ZD[IPI]-Z0 IC2MN=IC1+1 for (IC2=IC2MN; IC2<=NCP0; IC2++) { JIPC=JIPC0+IC2 IPI=IPC[JIPC] DX2=XD[IPI]-X0 DY2=YD[IPI]-Y0 DNMZ=DX1*DY2-DY1*DX2 if(DNMZ == 0.0) goto L22 DZ2=ZD[IPI]-Z0 DNMX=DY1*DZ2-DZ1*DY2 DNMY=DZ1*DX2-DX1*DZ2 if(DNMZ >= 0.0) goto L21 DNMX=-DNMX DNMY=-DNMY DNMZ=-DNMZ L21: NMX=NMX+DNMX NMY=NMY+DNMY NMZ=NMZ+DNMZ L22: } } JPD0=5*IP0 PD[JPD0-4]=-NMX/NMZ PD[JPD0-3]=-NMY/NMZ } //C ESTIMATION OF ZXX, ZXY, AND ZYY L30: for (IP0=1; IP0<=NDP0; IP0++) { JPD0=JPD0+5 X0=XD[IP0] JPD0=5*IP0 Y0=YD[IP0] ZX0=PD[JPD0-4] ZY0=PD[JPD0-3] NMXX=0.0 NMXY=0.0 NMYX=0.0 NMYY=0.0 NMZ =0.0 JIPC0=NCP0*(IP0-1) for (IC1=1; IC1<=NCPM1; IC1++) { JIPC=JIPC0+IC1 IPI=IPC[JIPC] DX1=XD[IPI]-X0 DY1=YD[IPI]-Y0 JPD=5*IPI DZX1=PD[JPD-4]-ZX0 DZY1=PD[JPD-3]-ZY0 IC2MN=IC1+1 for ( IC2=IC2MN; IC2<=NCP0; IC2++) { JIPC=JIPC0+IC2 IPI=IPC[JIPC] DX2=XD[IPI]-X0 DY2=YD[IPI]-Y0 DNMZ =DX1*DY2 -DY1*DX2 if(DNMZ == 0.0) goto L32 JPD=5*IPI DZX2=PD[JPD-4]-ZX0 DZY2=PD[JPD-3]-ZY0 DNMXX=DY1*DZX2-DZX1*DY2 DNMXY=DZX1*DX2-DX1*DZX2 DNMYX=DY1*DZY2-DZY1*DY2 DNMYY=DZY1*DX2-DX1*DZY2 if(DNMZ >= 0.0) goto L31 DNMXX=-DNMXX DNMXY=-DNMXY DNMYX=-DNMYX DNMYY=-DNMYY DNMZ =-DNMZ L31: NMXX=NMXX+DNMXX NMXY=NMXY+DNMXY NMYX=NMYX+DNMYX NMYY=NMYY+DNMYY NMZ =NMZ +DNMZ L32: } } PD[JPD0-2]=-NMXX/NMZ PD[JPD0-1]=-(NMXY+NMYX)/(2.0*NMZ) PD[JPD0] =-NMYY/NMZ } return } void idptip ( real scalar ndp, real vector xd, real vector yd, real vector zd, real scalar nt, real vector ipt, real scalar nl, real vector ipl, real vector pdd, real scalar iti, real scalar xii, real scalar yii, real scalar zii ) { /* ! !******************************************************************************* ! !! IDPTIP performs interpolation, determining a value of Z given X and Y. ! ! ! Modified: ! ! 19 February 2001 ! ! Parameters: ! ! Input, integer NDP, the number of data values. ! ! Input, real scalar XD(NDP), YD(NDP), the X and Y coordinates of the data. ! ! Input, real scalar ZD(NDP), the data values. ! ! Input, integer NT, the number of triangles. ! ! Input, ipt = integer array of dimension 3*nt containing the ! point numbers of the vertexes of the triangles, ! ! Input, integer NL, the number of border line segments. ! ! Input, integer IPL(3*NL), the point numbers of the end points of the ! border line segments and their respective triangle numbers, ! ! Input, real PDD(5*NDP). the partial derivatives at the data points, ! ! Input, integer ITI, triangle number of the triangle in which lies ! the point for which interpolation is to be performed, ! ! Input, real XII, YII, the X and Y coordinates of the point for which ! interpolation is to be performed. ! ! Output, real ZII, the interpolated Z value. ! */ // implicit none //real scalar ndp //real scalar nl //real scalar nt real scalar a real scalar aa real scalar ab real scalar ac real scalar act2 real scalar ad real scalar adbc //real scalar ap real scalar b real scalar bb real scalar bc real scalar bdt2 //real scalar bp real scalar c real scalar cc real scalar cd //real scalar cp real scalar csuv real scalar d real scalar dd real scalar dlt //real scalar dp real scalar dx real scalar dy real scalar g1 real scalar g2 real scalar h1 real scalar h2 real scalar h3 real scalar i real scalar idp real scalar il1 real scalar il2 //real scalar ipl(3*nl) //real scalar ipt(3*nt) real scalar it0 //real scalar iti //real scalar itpv real scalar jipl real scalar jipt real scalar jpd real scalar jpdd real scalar kpd real scalar ntl real scalar lu real scalar lv real scalar p0 //real scalar p00 //real scalar p01 //real scalar p02 //real scalar p03 //real scalar p04 //real scalar p05 real scalar p1 //real scalar p10 //real scalar p11 //real scalar p12 //real scalar p13 //real scalar p14 real scalar p2 //real scalar p20 //real scalar p21 //real scalar p22 //real scalar p23 real scalar p3 //real scalar p30 //real scalar p31 //real scalar p32 real scalar p4 //real scalar p40 //real scalar p41 //real scalar p5 //real scalar p50 real scalar pd pd=J(1,15,.) //real scalar pdd(5*ndp) real scalar thsv real scalar thus real scalar thuv real scalar thxu real scalar u real scalar v real scalar x x=J(1,3,.) //real scalar x0 //real scalar xd(*) //real scalar xii real scalar y y=J(1,3,.) //real scalar y0 //real scalar yd(*) //real scalar yii real scalar z z=J(1,3,.) real scalar z0 //real scalar zd(*) //real scalar zii real scalar zu zu=J(1,3,.) real scalar zuu zuu=J(1,3,.) real scalar zuv zuv=J(1,3,.) real scalar zv zv=J(1,3,.) real scalar zvv zvv=J(1,3,.) // save /idpt/ /* common /idpt/ itpv,x0,y0,ap,bp,cp,dp, & p00,p10,p20,p30,p40,p50,p01,p11,p21,p31,p41, & p02,p12,p22,p32,p03,p13,p23,p04,p14,p05 */ external real scalar itpv,x0,y0,ap,bp,cp,dp, p00,p10,p20,p30,p40,p50,p01,p11,p21,p31,p41, p02,p12,p22,p32,p03,p13,p23,p04,p14,p05 // Preliminary processing it0 = iti ntl = nt+nl if ( it0 > ntl ) { il1 = trunc(it0/ntl) il2 = it0-il1*ntl if(il1==il2) goto L40 goto L60 } // Calculation of ZII by interpolation. // Check if the necessary coefficients have been calculated. if ( it0 == itpv ) goto L30 // Load coordinate and partial derivative values at the vertexes. jipt = 3*(it0-1) jpd = 0 for (i = 1; i<=3; i++) { jipt = jipt+1 idp = ipt[jipt] x[i] = xd[idp] y[i] = yd[idp] z[i] = zd[idp] jpdd = 5*(idp-1) for (kpd = 1; kpd<=5; kpd++) { jpd = jpd+1 jpdd = jpdd+1 pd[jpd] = pdd[jpdd] } } // Determine the coefficients for the coordinate system // transformation from the XY system to the UV system and vice versa. x0 = x[1] y0 = y[1] a = x[2]-x0 b = x[3]-x0 c = y[2]-y0 d = y[3]-y0 ad = a*d bc = b*c dlt = ad-bc ap = d/dlt bp = -b/dlt cp = -c/dlt dp = a/dlt // Convert the partial derivatives at the vertexes of the // triangle for the UV coordinate system. aa = a*a act2 = 2.0E+00 *a*c cc = c*c ab = a*b adbc = ad+bc cd = c*d bb = b*b bdt2 = 2.0E+00 *b*d dd = d*d for (i = 1; i<=3; i++) { jpd = 5*i zu[i] = a*pd[jpd-4]+c*pd[jpd-3] zv[i] = b*pd[jpd-4]+d*pd[jpd-3] zuu[i] = aa*pd[jpd-2]+act2*pd[jpd-1]+cc*pd[jpd] zuv[i] = ab*pd[jpd-2]+adbc*pd[jpd-1]+cd*pd[jpd] zvv[i] = bb*pd[jpd-2]+bdt2*pd[jpd-1]+dd*pd[jpd] } // Calculate the coefficients of the polynomial. p00 = z[1] p10 = zu[1] p01 = zv[1] p20 = 0.5E+00 * zuu[1] p11 = zuv[1] p02 = 0.5E+00 * zvv[1] h1 = z[2]-p00-p10-p20 h2 = zu[2]-p10-zuu[1] h3 = zuu[2]-zuu[1] p30 = 10.0E+00 * h1 - 4.0E+00 * h2 + 0.5E+00 * h3 p40 = -15.0E+00 * h1 + 7.0E+00 * h2 - h3 p50 = 6.0E+00 * h1 - 3.0E+00 * h2 + 0.5E+00 * h3 h1 = z[3]-p00-p01-p02 h2 = zv[3]-p01-zvv[1] h3 = zvv[3]-zvv[1] p03 = 10.0E+00 * h1 - 4.0E+00 * h2 + 0.5E+00 * h3 p04 = -15.0E+00 * h1 + 7.0E+00 * h2 -h3 p05 = 6.0E+00 * h1 - 3.0E+00 * h2 + 0.5E+00 * h3 lu = sqrt(aa+cc) lv = sqrt(bb+dd) // thxu = atan2(c,a) thxu = atan2(a,c) // thuv = atan2(d,b)-thxu thuv = atan2(b,d)-thxu csuv = cos(thuv) p41 = 5.0E+00*lv*csuv/lu*p50 p14 = 5.0E+00*lu*csuv/lv*p05 h1 = zv[2]-p01-p11-p41 h2 = zuv[2]-p11-4.0E+00 * p41 p21 = 3.0E+00 * h1-h2 p31 = -2.0E+00 * h1+h2 h1 = zu[3]-p10-p11-p14 h2 = zuv[3]-p11- 4.0E+00 * p14 p12 = 3.0E+00 * h1-h2 p13 = -2.0E+00 * h1+h2 thus = atan2(b-a,d-c)-thxu thsv = thuv-thus aa = sin(thsv)/lu bb = -cos(thsv)/lu cc = sin(thus)/lv dd = cos(thus)/lv ac = aa*cc ad = aa*dd bc = bb*cc g1 = aa * ac*(3.0E+00*bc+2.0E+00*ad) g2 = cc * ac*(3.0E+00*ad+2.0E+00*bc) h1 = -aa*aa*aa*(5.0E+00*aa*bb*p50+(4.0E+00*bc+ad)*p41)-cc*cc*cc*(5.0E+00*cc*dd*p05+(4.0E+00*ad+bc)*p14) h2 = 0.5E+00 * zvv[2]-p02-p12 h3 = 0.5E+00 * zuu[3]-p20-p21 p22 = (g1*h2+g2*h3-h1)/(g1+g2) p32 = h2-p22 p23 = h3-p22 itpv = it0 // Convert XII and YII to UV system. L30: dx = xii-x0 dy = yii-y0 u = ap*dx+bp*dy v = cp*dx+dp*dy // Evaluate the polynomial. p0 = p00+v*(p01+v*(p02+v*(p03+v*(p04+v*p05)))) p1 = p10+v*(p11+v*(p12+v*(p13+v*p14))) p2 = p20+v*(p21+v*(p22+v*p23)) p3 = p30+v*(p31+v*p32) p4 = p40+v*p41 p5 = p50 zii = p0+u*(p1+u*(p2+u*(p3+u*(p4+u*p5)))) return // Calculation of ZII by extrapolation in the rectangle. // Check if the necessary coefficients have been calculated. L40: if(it0==itpv) goto L50 // Load coordinate and partial derivative values at the end // points of the border line segment. jipl = 3*(il1-1) jpd = 0 for (i = 1; i<=2; i++) { jipl = jipl+1 idp = ipl[jipl] x[i] = xd[idp] y[i] = yd[idp] z[i] = zd[idp] jpdd = 5*(idp-1) for (kpd = 1; kpd<=5; kpd++) { jpd = jpd+1 jpdd = jpdd+1 pd[jpd] = pdd[jpdd] } } // Determine the coefficients for the coordinate system // transformation from the XY system to the UV system // and vice versa. x0 = x[1] y0 = y[1] a = y[2]-y[1] b = x[2]-x[1] c = -b d = a ad = a*d bc = b*c dlt = ad-bc ap = d/dlt bp = -b/dlt cp = -bp dp = ap // Convert the partial derivatives at the end points of the // border line segment for the UV coordinate system. aa = a*a act2 = 2.0E+00 * a * c cc = c*c ab = a*b adbc = ad+bc cd = c*d bb = b*b bdt2 = 2.0E+00 * b * d dd = d*d for (i = 1; i<=2; i++) { jpd = 5*i zu[i] = a*pd[jpd-4]+c*pd[jpd-3] zv[i] = b*pd[jpd-4]+d*pd[jpd-3] zuu[i] = aa*pd[jpd-2]+act2*pd[jpd-1]+cc*pd[jpd] zuv[i] = ab*pd[jpd-2]+adbc*pd[jpd-1]+cd*pd[jpd] zvv[i] = bb*pd[jpd-2]+bdt2*pd[jpd-1]+dd*pd[jpd] } // Calculate the coefficients of the polynomial. p00 = z[1] p10 = zu[1] p01 = zv[1] p20 = 0.5E+00 * zuu[1] p11 = zuv[1] p02 = 0.5E+00 * zvv[1] h1 = z[2]-p00-p01-p02 h2 = zv[2]-p01-zvv[1] h3 = zvv[2]-zvv[1] p03 = 10.0E+00 * h1 - 4.0E+00*h2+0.5E+00*h3 p04 = -15.0E+00 * h1 + 7.0E+00*h2 -h3 p05 = 6.0E+00 * h1 - 3.0E+00*h2+0.5E+00*h3 h1 = zu[2]-p10-p11 h2 = zuv[2]-p11 p12 = 3.0E+00*h1-h2 p13 = -2.0E+00*h1+h2 p21 = 0.0E+00 p23 = -zuu[2]+zuu[1] p22 = -1.5E+00*p23 itpv = it0 // Convert XII and YII to UV system. L50: dx = xii-x0 dy = yii-y0 u = ap*dx+bp*dy v = cp*dx+dp*dy // Evaluate the polynomial. p0 = p00+v*(p01+v*(p02+v*(p03+v*(p04+v*p05)))) p1 = p10+v*(p11+v*(p12+v*p13)) p2 = p20+v*(p21+v*(p22+v*p23)) zii = p0+u*(p1+u*p2) return // Calculation of ZII by extrapolation in the triangle. // Check if the necessary coefficients have been calculated. L60: if ( it0 != itpv ) { // Load coordinate and partial derivative values at the vertex of the triangle. jipl = 3*il2-2 idp = ipl[jipl] x0 = xd[idp] y0 = yd[idp] z0 = zd[idp] jpdd = 5*(idp-1) for (kpd = 1; kpd<=5; kpd++) { jpdd = jpdd+1 pd[kpd] = pdd[jpdd] } // Calculate the coefficients of the polynomial. p00 = z0 p10 = pd[1] p01 = pd[2] p20 = 0.5E+00*pd[3] p11 = pd[4] p02 = 0.5E+00*pd[5] itpv = it0 } // Convert XII and YII to UV system. u = xii-x0 v = yii-y0 // Evaluate the polynomial. p0 = p00+v*(p01+v*p02) p1 = p10+v*p11 zii = p0+u*(p1+u*p20) return } void idsfft ( real scalar md, real scalar ncp, real scalar ndp, real vector xd, real vector yd, real vector zd, real scalar nxi, real scalar nyi, real scalar nzi, real vector xi, real vector yi, real matrix zi, real vector iwk, real vector wk ) { /* ! !******************************************************************************* ! !! IDSFFT fits a smooth surface Z(X,Y) given irregular (X,Y,Z) data. ! ! ! Discussion: ! ! IDSFFT performs smooth surface fitting when the projections of the ! data points in the (X,Y) plane are irregularly distributed. ! ! Special conditions: ! ! Inadequate work space IWK and WK may may cause incorrect results. ! ! The data points must be distinct and their projections in the XY ! plane must not be collinear, otherwise an error return occurs. ! ! Parameters: ! ! Input, integer MD, mode of computation (must be 1, 2, or 3, ! else an error return will occur). ! ! 1, if this is the first call to this routine, or if the value of ! NDP has been changed from the previous call, or if the contents of ! the XD or YD arrays have been changed from the previous call. ! ! 2, if the values of NDP and the XD, YD arrays are unchanged from ! the previous call, but new values for XI, YI are being used. If ! MD = 2 and NDP has been changed since the previous call to IDSFFT, ! an error return occurs. ! ! 3, if the values of NDP, NXI, NYI, XD, YD, XI, YI are unchanged ! from the previous call, i.e. if the only change on input to idsfft ! is in the ZD array. If MD = 3 and NDP, nxi or nyi has been changed ! since the previous call to idsfft, an error return occurs. ! ! Between the call with MD = 2 or MD = 3 and the preceding call, the ! iwk and wk work arrays should not be disturbed. ! ! Input, integer NDP, the number of data points. NDP must be at least 4. ! ! Input, real XD(NDP), YD(NDP), the X and Y coordinates of the data. ! ! Input, real ZD(NDP), the data values. ! ! Input, integer NXI, NYI, the number of output grid points in the ! X and Y directions. NXI and NYI must each be at least 1. ! ! Input, integer NZI, the first dimension of ZI. NZI must be at ! least NXI. ! ! Input, real XI(NXI), YI(NYI), the X and Y coordinates of the grid ! points. ! ! Workspace, integer IWK(31*NDP+NXI*NYI). ! ! Workspace, real WK(6*NDP). ! ! Output, real ZI(NZI,NYI), contains the interpolated Z values at the ! grid points. ! */ // implicit none //real scalar ndp //real scalar nxi //real scalar nyi //real scalar nzi //real scalar ap //real scalar bp //real scalar cp //real scalar dp real scalar il1 real scalar il2 real scalar iti //real scalar itpv //real scalar iwk(31*ndp + nxi*nyi) real scalar ixi real scalar iyi real scalar izi real scalar jig0mn real scalar jig0mx real scalar jig1mn real scalar jig1mx real scalar jigp real scalar jngp real scalar jwigp real scalar jwigp0 real scalar jwipl real scalar jwipt real scalar jwiwl real scalar jwiwp real scalar jwipc real scalar jwngp real scalar jwngp0 real scalar jwwpd //real scalar md real scalar ngp0 real scalar ngp1 real scalar nl real scalar nngp real scalar nt //real scalar p00 //real scalar p01 //real scalar p02 //real scalar p03 //real scalar p04 //real scalar p05 //real scalar p10 //real scalar p11 //real scalar p12 //real scalar p13 //real scalar p14 //real scalar p20 //real scalar p21 //real scalar p22 //real scalar p23 //real scalar p30 //real scalar p31 //real scalar p32 //real scalar p40 //real scalar p41 //real scalar p50 //real scalar wk(6*ndp) //real scalar x0 //real scalar xd(ndp) //real scalar xi(nxi) //real scalar y0 //real scalar yd(ndp) //real scalar yi(nyi) //real scalar zd(ndp) //real scalar zi(nzi,nyi) real vector ipt, ipl, iwl, iwp, ipc real vector ngp, igp real vector wpd real scalar tzi /* save /idpt/ ! common /idpt/ itpv,x0,y0,ap,bp,cp,dp, & p00,p10,p20,p30,p40,p50,p01,p11,p21,p31,p41, & p02,p12,p22,p32,p03,p13,p23,p04,p14,p05 */ external real scalar itpv,x0,y0,ap,bp,cp,dp, p00,p10,p20,p30,p40,p50,p01,p11,p21,p31,p41, p02,p12,p22,p32,p03,p13,p23,p04,p14,p05 // Error check. if ( md < 1 | md > 3 ) { _error(6,"IDSFFT - Fatal error: Input parameter MD out of range.") return } if (ncp!=-1 & (ncp < 2 | ncp>min((25,ndp))) ) { _error(10,"IDSFFT - Fatal error: Input parameter NCP out of range.") return } if ( ndp < 4 ) { _error(7,"IDSFFT - Fatal error: Input parameter NDP out of range.") return } if ( nxi < 1 | nyi < 1 ) { _error(8,"IDSFFT - Fatal error: Input parameter NXI or NYI out of range.") return } if ( nxi > nzi ) { _error(9,"IDSFFT - Fatal error: Input parameter NZI is less than NXI.") return } if ( md <= 1 ) { iwk[1] = ncp iwk[2] = ndp } else { if ( ncp != iwk[1] ) { _error(10,"IDSFFT - Fatal error: MD = 2 or 3 but ncp was changed since last call.") return } if ( ndp != iwk[2] ) { _error(11,"IDSFFT - Fatal error: MD = 2 or 3 but ndp was changed since last call.") return } } if ( md <= 2 ) { iwk[3] = nxi iwk[4] = nyi } else { if ( nxi != iwk[3] ) { _error(12,"IDSFFT - Fatal error: MD = 3 but nxi was changed since last call.") return } if ( nyi != iwk[4] ) { _error(13,"IDSFFT - Fatal error: MD = 3 but nyi was changed since last call.") return } } // Allocation of storage areas in the IWK array. jwipt = 16 jwiwl = 6*ndp+1 jwngp0 = jwiwl-1 jwipl = 24*ndp+1 jwiwp = 30*ndp+1 //jwigp0 = 31*ndp jwwpd = 5*ndp+1 jwipc=27*ndp+1 jwigp0=max((31,27+ncp))*ndp // Triangulate the XY plane. if ( md == 1 ) { ipt=iwk[|jwipt\.|] ipl=iwk[|jwipl\.|] iwl=iwk[|jwiwl\.|] iwp=iwk[|jwiwp\.|] // idtang ( ndp, xd, yd, nt, iwk[|jwipt\.|], nl, iwk[|jwipl\.|], iwk[|jwiwl\.|], iwk[|jwiwp\.|], wk ) idtang ( ndp, xd, yd, nt, ipt, nl, ipl, iwl, iwp, wk ) iwk[|jwipt\.|]=ipt iwk[|jwiwl\.|]=iwl // Why not done before? iwk[|jwipl\.|]=ipl iwk[|jwiwp\.|]=iwp // Why not done before? iwk[5] = nt iwk[6] = nl if ( nt == 0 ) { return } } else { nt = iwk[5] nl = iwk[6] } // DETERMINES NCP POINTS CLOSEST TO EACH DATA POINT. (FOR MD=1) // 50 IF(MD0.GT.1) GO TO 60 // CALL IDCLDP(NDP0,XD,YD,NCP0,IWK(JWIPC)) // IF(IWK(JWIPC).EQ.0) RETURN if (md == 1 & ncp>=2) { idcldp(ndp, xd, yd, ncp, ipc) //iwk[|jwipc\.|]=ipc // For some reason this doesn't work here for (i=1; i<=cols(ipc); i++) { iwk[jwipc+i-1]=ipc[i] } if (iwk[jwipc]==0) { return } } // Sort output grid points in ascending order of the triangle // number and the border line segment number. if ( md <= 2 ) { ipt=iwk[|jwipt\.|] ipl=iwk[|jwipl\.|] ngp=iwk[|jwngp0+1\.|] igp=iwk[|jwigp0+1\.|] // idgrid ( xd, yd, nt, iwk[|jwipt\.|], nl, iwk[|jwipl\.|], nxi, nyi, xi, yi, iwk[|jwngp0+1\.|], iwk[|jwigp0+1\.|] ) idgrid ( xd, yd, nt, ipt, nl, ipl, nxi, nyi, xi, yi, ngp, igp ) iwk[|jwngp0+1\.|]=ngp iwk[|jwigp0+1\.|]=igp } // Estimate partial derivatives at all data points. // idpdrv ( ndp, xd, yd, zd, nt, iwk[|jwipt\.|], wk, wk[|jwwpd\.|] ) if (ncp==-1) { wpd=wk[|jwwpd\.|] idpdrv ( ndp, xd, yd, zd, nt, ipt, wk, wpd ) wk[|jwwpd\.|]=wpd } else { ipc=iwk[|jwipc\.|] idpdrv_near ( ndp, xd, yd, zd, ncp, ipc, wk) iwk[|jwipc\.|]=ipc } // Interpolate the ZI values. itpv = 0 jig0mx = 0 jig1mn = nxi*nyi+1 nngp = nt+2*nl for (jngp = 1; jngp<=nngp; jngp++) { iti = jngp if ( jngp > nt ) { il1 = trunc((jngp-nt+1)/2) il2 = trunc((jngp-nt+2)/2) if(il2>nl) { il2 = 1 } iti = il1*(nt+nl)+il2 } jwngp = jwngp0+jngp ngp0 = iwk[jwngp] if ( ngp0 != 0 ) { jig0mn = jig0mx+1 jig0mx = jig0mx+ngp0 for (jigp = jig0mn; jigp<=jig0mx; jigp++) { jwigp = jwigp0+jigp izi = iwk[jwigp] iyi = trunc((izi-1)/nxi)+1 ixi = izi-nxi*(iyi-1) ipt=iwk[|jwipt\.|] ipl=iwk[|jwipl\.|] tzi=. // idptip ( ndp, xd, yd, zd, nt, iwk[|jwipt\.|], nl, iwk[|jwipl\.|], wk, iti, xi[ixi], yi[iyi], zi[ixi,iyi] ) idptip ( ndp, xd, yd, zd, nt, ipt, nl, ipl, wk, iti, xi[ixi], yi[iyi], tzi ) zi[ixi,iyi]=tzi } } jwngp = jwngp0+2*nngp+1-jngp ngp1 = iwk[jwngp] if ( ngp1 != 0 ) { jig1mx = jig1mn-1 jig1mn = jig1mn-ngp1 for (jigp = jig1mn; jigp<=jig1mx; jigp++) { jwigp = jwigp0+jigp izi = iwk[jwigp] iyi = trunc((izi-1)/nxi+1) ixi = izi-nxi*(iyi-1) ipt=iwk[|jwipt\.|] ipl=iwk[|jwipl\.|] tzi=. // idptip ( ndp, xd, yd, zd, nt, iwk[|jwipt\.|], nl, iwk[|jwipl\.|], wk, iti, xi[ixi], yi[iyi], zi[ixi,iyi] ) idptip ( ndp, xd, yd, zd, nt, ipt, nl, ipl, wk, iti, xi[ixi], yi[iyi], tzi ) zi[ixi,iyi]=tzi } } } return } // Second version of spdt() that is used only by IDTANG real scalar t_spdt(real scalar u1, real scalar v1, real scalar u2, real scalar v2,real scalar u3,real scalar v3) { return((u2-u1)*(u3-u1)+(v2-v1)*(v3-v1)) } // Second version of vpdt() that is used only by IDTANG real scalar t_vpdt(real scalar u1, real scalar v1, real scalar u2, real scalar v2,real scalar u3,real scalar v3) { return((v3-v1)*(u2-u1)-(u3-u1)*(v2-v1)) } void idtang ( real scalar ndp, real vector xd, real vector yd, real scalar nt, real vector ipt, real scalar nl, real vector ipl, real vector iwl, real vector iwp, real vector wk ) { /* ! !******************************************************************************* ! !! IDTANG performs triangulation. ! ! ! Discussion: ! ! The routine divides the XY plane into a number of triangles according to ! given data points in the plane, determines line segments that form ! the border of data area, and determines the triangle numbers ! corresponding to the border line segments. ! ! At completion, point numbers of the vertexes of each triangle ! are listed counter-clockwise. Point numbers of the end points ! of each border line segment are listed counter-clockwise, ! listing order of the line segments being counter-clockwise. ! ! Parameters: ! ! Input, integer NDP, the number of data points. ! ! Input, real XD(NDP), YD(NDP), the X and Y coordinates of the data. ! ! Output, integer NT, the number of triangles, ! ! Output, integer IPT(6*NDP-15), where the point numbers of the ! vertexes of the IT-th triangle are to be stored as entries ! 3*IT-2, 3*IT-1, and 3*IT, for IT = 1 to NT. ! ! Output, integer NL, the number of border line segments. ! ! Output, integer IPL(6*NDP), where the point numbers of the end ! points of the (il)th border line segment and its respective triangle ! number are to be stored as the (3*il-2)nd, (3*il-1)st, and (3*il)th ! elements, il = 1,2,..., nl. ! ! Workspace, integer IWL(18*NDP), ! ! Workspace, integer IWP(NDP), ! ! Workspace, real WK(NDP). ! */ // implicit none //integer ndp //real scalar dsqf real scalar dsqi real scalar dsqmn real scalar epsln epsln = 1.0E-06 //real scalar idxchg real scalar il real scalar ilf real scalar iliv real scalar ilt3 real scalar ilvs real scalar ip real scalar ip1 real scalar ip1p1 real scalar ip2 real scalar ip3 //real scalar ipl(6*ndp) real scalar ipl1 real scalar ipl2 real scalar iplj1 real scalar iplj2 real scalar ipmn1 real scalar ipmn2 //real scalar ipt(6*ndp-15) real scalar ipt1 real scalar ipt2 real scalar ipt3 real scalar ipti real scalar ipti1 real scalar ipti2 real scalar irep real scalar it real scalar it1t3 real scalar it2t3 real vector itf itf=J(1,2,.) real scalar its real scalar itt3 real scalar itt3r //real scalar iwl(18*ndp) //real scalar iwp(ndp) real scalar ixvs real scalar ixvspv real scalar jl1 real scalar jl2 real scalar jlt3 real scalar jp real scalar jp1 real scalar jp2 real scalar jpc real scalar jpmn real scalar jpmx real scalar jwl real scalar jwl1 real scalar jwl1mn //real scalar nl real scalar nl0 real scalar nlf real scalar nlfc real scalar nlft2 real scalar nln real scalar nlnt3 real scalar nlsh real scalar nlsht3 real scalar nlt3 real scalar nrep nrep = 100 //real scalar nt real scalar nt0 real scalar ntf real scalar ntt3 real scalar ntt3p3 real scalar sp //real scalar spdt //real scalar u1 //real scalar u2 //real scalar u3 //real scalar v1 //real scalar v2 //real scalar v3 real scalar vp //real scalar vpdt //real scalar wk(ndp) real scalar x1 real scalar x2 real scalar x3 //real scalar xd(ndp) real scalar xdmp real scalar y1 real scalar y2 real scalar y3 //real scalar yd(ndp) real scalar ydmp /* ! ! Statement functions ! dsqf(u1,v1,u2,v2) = (u2-u1)**2+(v2-v1)**2 spdt(u1,v1,u2,v2,u3,v3) = (u2-u1)*(u3-u1)+(v2-v1)*(v3-v1) vpdt(u1,v1,u2,v2,u3,v3) = (v3-v1)*(u2-u1)-(u3-u1)*(v2-v1) */ // Preliminary processing if ( ndp < 4 ) { _error(12,"IDTANG - Fatal error: Input parameter NDP out of range.") return } // Determine IPMN1 and IPMN2, the closest pair of data points. dsqmn = dsqf(xd[1],yd[1],xd[2],yd[2]) ipmn1 = 1 ipmn2 = 2 for (ip1 = 1; ip1<=ndp-1; ip1++) { x1 = xd[ip1] y1 = yd[ip1] ip1p1 = ip1+1 for (ip2 = ip1p1; ip2<=ndp; ip2++) { dsqi = dsqf(x1,y1,xd[ip2],yd[ip2]) if ( dsqi == 0.0 ) { _error(13,"IDTANG - Fatal error: Two of the input data points are identical.") return } if(dsqi ( abs(sp) * epsln ) ) goto L37 } _error(14,"IDTANG - Fatal error: All collinear data points.") return L37: if ( jp != 3 ) { jpmx = jp for (jpc = 4; jpc<=jpmx; jpc++) { jp = jpmx+4-jpc iwp[jp] = iwp[jp-1] } iwp[3] = ip } // Form the first triangle. // Store point numbers of the vertexes of the triangle in the IPT array, // store point numbers of the border line segments and the triangle number in // the IPL array. ip1 = ipmn1 ip2 = ipmn2 ip3 = iwp[3] if ( t_vpdt(xd[ip1],yd[ip1],xd[ip2],yd[ip2],xd[ip3],yd[ip3]) < 0.0E+00 ) { ip1 = ipmn2 ip2 = ipmn1 } nt0 = 1 ntt3 = 3 ipt[1] = ip1 ipt[2] = ip2 ipt[3] = ip3 nl0 = 3 nlt3 = 9 ipl[1] = ip1 ipl[2] = ip2 ipl[3] = 1 ipl[4] = ip2 ipl[5] = ip3 ipl[6] = 1 ipl[7] = ip3 ipl[8] = ip1 ipl[9] = 1 // Add the remaining data points, one by one. for (jp1 = 4; jp1<=ndp; jp1++) { ip1 = iwp[jp1] x1 = xd[ip1] y1 = yd[ip1] // Determine the first invisible and visible border line segments, iliv and // ilvs. for (il = 1; il <= nl0; il++) { ip2 = ipl[3*il-2] ip3 = ipl[3*il-1] x2 = xd[ip2] y2 = yd[ip2] x3 = xd[ip3] y3 = yd[ip3] sp = t_spdt(x1,y1,x2,y2,x3,y3) vp = t_vpdt(x1,y1,x2,y2,x3,y3) if ( il == 1 ) { ixvs = 0 if(vp<=(abs(sp)*(-epsln))) ixvs = 1 iliv = 1 ilvs = 1 goto L53 } ixvspv = ixvs if ( vp <= (abs(sp)*(-epsln)) ) { ixvs = 1 if(ixvspv==1) goto L53 ilvs = il if(iliv!=1) goto L54 goto L53 } ixvs = 0 if ( ixvspv != 0 ) { iliv = il if(ilvs!=1) goto L54 } L53: } if(iliv==1 & ilvs==1) ilvs = nl0 L54: if(ilvs 0.0E+00 ) { u1 = (y3-y1)*(x4-x1)-(x3-x1)*(y4-y1) u2 = (y4-y2)*(x3-x2)-(x4-x2)*(y3-y2) a1sq = (x1-x3)^2+(y1-y3)^2 a4sq = (x4-x1)^2+(y4-y1)^2 c1sq = (x3-x4)^2+(y3-y4)^2 a2sq = (x2-x4)^2+(y2-y4)^2 a3sq = (x3-x2)^2+(y3-y2)^2 c3sq = (x2-x1)^2+(y2-y1)^2 s1sq = u1*u1 / (c1sq*max((a1sq,a4sq))) s2sq = u2*u2 / (c1sq*max((a2sq,a3sq))) s3sq = u3*u3 / (c3sq*max((a3sq,a1sq))) s4sq = u4*u4 / (c3sq*max((a4sq,a2sq))) if ( min (( s3sq, s4sq )) - min (( s1sq, s2sq )) > epsln ) { idx = 1 } } idxchg = idx return (idxchg) } end // My version of cvxhull (Same as the original but just included here so it doesn't have to be bundled) * 1.1 RAR 21 January 2004, with improved coding * 1.0.1 RAR 12 December 2003, with suggestions from NJC * 1.0 RAR 22 November 2003 - England won RU World Cup program bip_cvxhull, sortpreserve version 8.0 syntax varlist(min=2 max=2) [if] [in] /// [ , Hulls(numlist int min=0 max=2 >0) MDPlot(int 8) /// GROup(varname) SELect(numlist int min=0 >0 sort) /// MEAns PREfix(string) noREPort noRETain /// SCATopt(string) noGRAph SAVing(string)] loc reporting = ("`report'" != "noreport") if "`retain'" == "noretain" & "`graph'" == "nograph" { di as err "That combination of options is silly: quitting now" exit 498 } * Nick Cox thinks this parsing is still flakey if `"`saving'"' != "" { // Check saving filename tokenize `"`saving'"', parse(",") args saving comma replace if index(`"`saving'"', ".") == 0 loc saving `"`saving'.gph"' capture confirm file `saving' if _rc == 0 { /* existing file */ if "`replace'" != "replace" { di as err "Invalid save file - did you mean to replace?" exit 198 } loc saving `"`saving',replace"' } else confirm new file `saving' // invalid file - unless new loc saving `"saving("`saving'")"' } if "`hulls'" == "" { loc hulls = 1 loc hullgap = 1 // Set default values } else { tokenize "`hulls'" loc hulls = `1' if "`2'" == "" loc hullgap = 1 else loc hullgap = `2' } marksample touse // Deal with `if', `in' & missing qui count if `touse' if r(N) == 0 error 2000 tempvar sample grp scratch count pts hullno qui gen `sample'=`touse' // & remember all points considered if "`group'" == "" { // group variable may be absent, string or numeric qui gen `grp' = 1 loc select "1" loc maxgrp = 1 } else { qui egen `grp' = group(`group') if `touse', label if "`select'" != "" { loc j = 0 foreach i in `select' { loc ++j } loc maxgrp = `j' // no of items } else { su `grp', meanonly loc maxgrp = r(max) loc select "1 / `maxgrp'" } qui { egen `scratch' = eqany(`grp'), values(`select') replace `touse' = `touse' * `scratch' } if `reporting' di as txt "Codes for groups based on " as res "`group'" label list `grp' } * Double check, in case group selection has dropped all obs qui count if `touse' if r(N) == 0 error 2000 loc retd = 1 + int((`hulls'-1)/`hullgap') // number of hulls loc hull = plural(`retd', "hull") loc gs = plural(`maxgrp', "group") if `reporting' di as txt "Up to `retd' `hull' to be saved for `maxgrp' `gs'" if "`prefix'" == "" loc prefix "_cvxh" capture drop `prefix'*l capture drop `prefix'*r capture drop `prefix'grp capture drop `prefix'hull capture drop `prefix'pts capture drop `prefix'cnt capture drop `prefix'*mindex capture drop `prefix'*maxdex capture drop `scratch' tokenize `varlist' args y x sort `grp' `x' `y' // within group, by ascending x then y qui { gen `count' = 0 gen `pts' = 0 gen `hullno' = 0 } tempname gap loc gap = `hullgap' -1 // force retain hull 1 loc maxhull = 0 ***** Main loop ****************************************** tempvar leftpath rightpath onhull // re-use scratch for points on segment qui { gen `leftpath' = . // set up outside loop to use replace gen `rightpath' = . // within loop & save RAM allocation time gen `onhull' = 0 gen `scratch' = 0 } forvalues h = 1 / `hulls' { // b1 : brackets numbered for sanity! qui count if `touse' // stop when all points peeled if r(N) > 0 { // b2 qui { replace `leftpath' = . replace `rightpath' = . replace `onhull' = 0 } loc sp = 1 // Starting point loc notstarted = 1 while `sp' <= _N { // b3 scan all observations loc curgrp = `grp'[`sp'] if !`touse'[`sp'] { loc ++sp // ignore point & increment pointer } else { // b4 - point processing ... while `grp'[`sp'] == `curgrp' { // b5 within group if `notstarted' { // b6 mark first point of current group qui { replace `leftpath' = `y' in `sp' replace `rightpath' = `y' in `sp' replace `onhull' = (_n==`sp') replace `scratch' = 0 } loc leftcur = `sp' loc rightcur = `sp' loc sp1 = `sp' loc curgrp = `grp'[`sp'] loc mindex = 1 loc notstarted = 0 loc ++sp } else { // 6 loc j = `leftcur' + 1 loc maxL = -1 loc maxD = 0 loc found = 0 while `j'<=_N & `grp'[`j']==`curgrp' { // b7 find next left */ if `touse'[`j']==1 { // b8 loc d = sqrt( (`y'[`j'] - `y'[`leftcur'])^2 + (`x'[`j'] - `x'[`leftcur'])^2 ) if float(`d') > 0 { // b9 loc cosa = (`y'[`j'] - `y'[`leftcur'])/ `d' if float(`cosa') > float(`maxL') { // b10 new leftmost direction loc maxL = `cosa' loc next = `j' qui replace `scratch' = (_n==`j') qui replace `scratch' = 1 in `sp1' loc found = 1 } else { if float(`cosa') == float(`maxL') { // collinear if float(`d') > float(`maxD') loc next = `j' qui replace `scratch' = 1 in `j' loc found = 1 } } } // 9 else { // 9 qui replace `onhull' = 1 in `j' // coincident with current point loc next = `j' } //9 } // 8 loc ++j } // 7 qui replace `onhull' = (`onhull' | `scratch') if `found' { qui { replace `leftpath' = `y' in `next' loc ++mindex loc leftcur = `next' } } loc j = `rightcur' + 1 loc maxR = 1 loc maxD = 0 loc found = 0 while `j'<=_N & `grp'[`j']==`curgrp' { // b7 find next right if `touse'[`j']==1 { //8 loc d = sqrt( (`y'[`j'] - `y'[`rightcur'])^2 + (`x'[`j'] - `x'[`rightcur'])^2 ) if float(`d') > 0 { //9 loc cosa = (`y'[`j'] - `y'[`rightcur']) / `d' if float(`cosa') < float(`maxR') { //10 loc maxR = `cosa' loc next = `j' qui replace `scratch' = (_n==`j') qui replace `scratch' = 1 in `sp1' loc found = 1 } else { if float(`cosa') == float(`maxR') { // collinear if float(`d') > float(`maxD') loc next = `j' qui replace `scratch' = 1 in `j' loc found = 1 } } } //9 else { //9 qui replace `onhull' = 1 in `j' // coincident loc next = `j' } //9 } //8 loc ++j } // 7 qui replace `onhull' = (`onhull' | `scratch') if `found' { qui { replace `rightpath' = `y' in `next' loc ++mindex loc rightcur = `next' } } * Check if found a hull, or continue loc sp = `leftcur' if `rightcur' < `leftcur' loc sp = `rightcur' if `leftcur' == `rightcur' { // hull closed qui { replace `hullno' = -`h' in `sp' replace `hullno' = `h' in `sp1' replace `count' = `mindex' - 1 in `sp1' // last point double counted qui count if `onhull' replace `count' = r(N) in `sp' replace `pts' = `h' if `onhull' replace `touse' = 0 if `onhull' } loc notstarted = 1 while `grp'[`sp'] == `curgrp' { loc ++sp } } } // e6 } // e5 found end of group if `mindex' == 1 { // special case - single point hull qui { qui count if `onhull' replace `count' = r(N) in `sp1' replace `hullno' = `h' in `sp1' loc notstarted = 1 } } } // e4 end point processing } // e3 end loop over data * Save hull if requested loc ++gap if `gap' == `hullgap' { loc ++maxhull qui { gen `prefix'`h'l = `leftpath' gen `prefix'`h'r = `rightpath' } if `maxhull' <= `mdplot' loc hulllist "`prefix'1l-`prefix'`h'r" loc gap = 0 // restart count if `reporting' di as txt " hull level `h' calculated and saved" } } // e2 do nothing if no data points } // e1 exit loop when all hulls marked ***** Optional plot if `"`graph'"' != "nograph" { if `reporting' di as txt "Graph will be plotted presently" * Build and execute graph command ... as a very long macro text! loc colours "black black blue blue dkorange dkorange magenta magenta emerald emerald khaki khaki cyan cyan red red" * Set up point plot and add line plots for each selected group * RAR's preference for default horizontal y labels loc gr `"scatter `y' `x' if `sample',yti("`y'")yla(,angle(0))`saving' `scatopt'"' if "`means'" == "means" { tempvar ymean xmean egen `ymean' = mean(`y') if `select', by(`grp') egen `xmean' = mean(`x') if `select', by(`grp') loc gr `"`gr'||scatter `ymean' `xmean',ms(T)xti("`x'")"' } foreach i of numlist `select' { loc gr `"`gr'||line `hulllist' `x' if `grp'==`i',clc(`colours')legend(off)"' } `gr' // & execute macro command } if "`retain'" != "noretain" { qui { if "`group'" != "" gen `prefix'grp = `grp' gen `prefix'hull = `hullno' gen `prefix'cnt = `count' gen `prefix'pts = `pts' } } else { capture drop `prefix'*l `prefix'*r } if `reporting' di as txt "cvxhull run" end // of cvxhull