program main C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C % C Copyright (C) 1996, The Board of Trustees of the Leland Stanford % C Junior University. All rights reserved. % C % C The programs in GSLIB are distributed in the hope that they will be % C useful, but WITHOUT ANY WARRANTY. No author or distributor accepts % C responsibility to anyone for the consequences of using them or for % C whether they serve any particular purpose or work at all, unless he % C says so in writing. Everyone is granted permission to copy, modify % C and redistribute the programs in GSLIB, but only under the condition % C that this notice and the above copyright notice remain intact. % C % C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c----------------------------------------------------------------------- c c Indicator Variograms from BiGaussian Distribution c ************************************************* c c or c ** c c Normal Scores Variogram (of a Baussian Distribution) c **************************************************** c c from Indicator Variogram c ************************ c c c This program returns: c c the values of the theoretical indicator semivariograms for various c thresholds given a bivariate Gaussian distribution with a specified c input normal scores semivariogram c c or c c the values of the theoretical normal scores semivariogram of a c bivariate Gaussian distribution that, after truncation at the 1-p quantile, c will yield a binary image having the desired (input) standardized indicator c semivariogram (1-p being the input cutoff/proportion) c c c c INPUT/OUTPUT PARAMETERS: c c imode input mode: model (1), variogram file (2) c infl input variogram file c pcut proportion / cutoff c icalc indicator-->nscores (1) nsores-->indicator (2) c outfl the output file for vargplt c ncut number of cutoffs c cut() the cutoffs (in cdf units) c ndir,nlag number of directions, number of lags c xoff,yoff,zoff for each direction: the specification c nst,c0 Normal scores variogram: nst, nugget c it,aa,cc type, a parameter, c parameter c ang1,ang2,...,anis2 anisotropy definition c c c c PROGRAM NOTES: c c 1. Setting the number of cutoffs to a -1 times the number of c cutoffs will cause the variograms to be standardized to a sill c of one. c c 2. Only 1 cutoff/proportion is allowed when calculating ry from ri. c However, various directional variograms can be processed. c c 3. When imode=2, i.e., nscore --> indicator, then it is assumed that c the input variogram file is a sample variogram file output from c gamv or gam programs. Not a vmodel output file. c c c NO EXTERNAL REFERENCES c c c Phaedon C. Kyriakidis May 1997 c----------------------------------------------------------------------- parameter (PI=3.14159265, DEG2RAD=PI/180.0, EPSLON = 1.0e-20, + MLAG=100, MDIR=4, MAXNST=4, MAXROT=4, VERSION=2.00, + NTERMS=79) implicit double precision (a,r,g) dimension a(NTERMS),rnew(MLAG*MDIR),r(MLAG*MDIR), + g(MLAG*MDIR),cnew(MLAG*MDIR),gnew(MLAG*MDIR) real h(MLAG*MDIR),maxcov, gnewst(MLAG*MDIR),xoff(MDIR), + yoff(MDIR),zoff(MDIR),hm(MLAG*MDIR),tm(MLAG*MDIR) real c0(1),aa(MAXNST),cc(MAXNST),ang1(MAXNST),ang2(MAXNST), + ang3(MAXNST),anis1(MAXNST),anis2(MAXNST), zc(100) real*8 rotmat(MAXROT,3,3),pcut,xp integer nst(1),it(MAXNST) character infl*40,outfl*40,str*80 logical testfl,corr,first,finish data lin/1/,lout/2/,corr/.false./, + finish/.false./,first/.true./ c c Note VERSION number before anything else: c write(*,9999) VERSION 9999 format(/' BIGAUS2 Version: ',f5.3/) c c Get the name of the parameter file - try the default name if no input: c write(*,*) 'Which parameter file do you want to use?' read (*,'(a20)') str(1:20) if(str(1:1).eq.' ') str(1:20) = 'bigaus2.par ' inquire(file=str(1:20),exist=testfl) if(.not.testfl) then write(*,*) 'ERROR - the parameter file does not exist,' write(*,*) ' check for the file and try again ' write(*,*) if(str(1:20).eq.'bigaus2.par ') then write(*,*) ' creating a blank parameter file' call makepar write(*,*) end if stop endif open(lin,file=str(1:20),status='old') c c Find Start of Parameters: c 1 read(lin,'(a4)',end=98) str(1:4) if(str(1:4).ne.'STAR') go to 1 c c Read Input Parameters: c read(lin,*,err=98) imode write(*,*) ' input mode = ',imode read(lin,'(a40)',err=98) infl call chknam(infl,40) write(*,*) ' input file = ',infl if(imode.eq.2) then inquire(file=infl,exist=testfl) if(.not.testfl) then write(*,*) 'ERROR - imode=2: variogram file is required,' write(*,*) ' check for file and try again ' stop endif endif read(lin,*,err=98) zcut write(*,*) ' threshold probability = ',zcut read(lin,*,err=98) icalc write(*,*) ' calculation mode = ',icalc read(lin,'(a40)',err=98) outfl call chknam(outfl,40) write(*,*) ' output file = ',outfl read(lin,*,err=98) ncut write(*,*) ' number of cutoffs = ',ncut if(ncut.lt.0) then ncut = -ncut corr = .true. endif read(lin,*,err=98) (zc(i),i=1,ncut) write(*,*) ' cutoffs = ',(zc(i),i=1,ncut) read(lin,*,err=98) ndir,nlag write(*,*) ' ndir,nlag = ',ndir,nlag do i=1,ndir read(lin,*,err=98) azm,dip,xlag write(*,*) ' azm, dip, lag = ',azm,dip,xlag xoff(i) = sin(DEG2RAD*azm)*xlag*cos(DEG2RAD*dip) yoff(i) = cos(DEG2RAD*azm)*xlag*cos(DEG2RAD*dip) zoff(i) = sin(DEG2RAD*dip)*xlag write(*,*) ' x,y,z offsets = ',xoff(i),yoff(i),zoff(i) end do read(lin,*,err=98) nst(1),c0(1) write(*,*) ' nst, c0 = ',nst(1),c0(1) if(nst(1).le.0) then write(*,9996) nst(1) 9996 format(' nst must be at least 1, it has been set to ',i4,/, + ' The c or a values can be set to zero') stop endif do i=1,nst(1) read(lin,*,err=98) it(i),cc(i),ang1(i),ang2(i),ang3(i) read(lin,*,err=98) aa(i),aa1,aa2 anis1(i) = aa1 / max(aa(i),EPSLON) anis2(i) = aa2 / max(aa(i),EPSLON) if(it(i).eq.4) then if(aa(i).lt.0.0) stop ' INVALID power variogram' if(aa(i).gt.2.0) stop ' INVALID power variogram' end if write(*,*) ' it,cc,ang[1,2,3] = ',it(i),cc(i), + ang1(i),ang2(i),ang3(i) write(*,*) ' a1 a2 a3 = ',aa(i),aa1,aa2 end do close(lin) c c Main if statement (input mode) c if(imode.eq.2) then pcut = dble(zcut) p1p = dble(pcut*(1.0d0-pcut)) c c Loop through the variogram file: c open(lin,file=infl,status='OLD') read(lin,'(a80)',end=31) str open(lout,file=outfl,status='UNKNOWN') write(lout,*) str(1:40) iv = 1 if(icalc.eq.2) then call coeffs(pcut,a,NTERMS) endif 11 np = 0 22 read(lin,'(a80)',end=31) str c c Kludge: c go to 41 31 finish=.true. str(1:1) = 'a' 41 continue c c Process this line of input from variogram file: c if(str(1:1).ne.' ') then if(np.eq.0) go to 11 iv = iv + 1 c c Append this variogram to the current file? c if(first.eq..true..and.finish.ne..true.) then write(lout,*) str(1:40) first = .false. endif if(finish) then close(lout) close(lin) go to 9997 endif go to 11 else np = np + 1 backspace lin read(lin,*) ii,h(np),g(np),nprs,hm(np),tm(np) r(np) = dble(1.0d0-g(np)) c c Sort out which correlation coefficient is needed c if(icalc.eq.1) then c c Normal scores input ---> standardized indicator output c rnew(np) = rhoi(r(np),pcut,NTERMS) gnew(np) = 1.0d0-rnew(np) write(lout,59) ii,h(np),gnew(np),nprs,hm(np),tm(np) else c c Standardized indicator input ---> normal scores output c rnew(np) = rhoy(r(np),pcut,a,NTERMS) gnew(np) = 1.0d0-rnew(np) write(lout,59) ii,h(np), gnew(np),nprs,hm(np),tm(np) endif 59 format(1x,i3,1x,f12.3,1x,f12.5,1x,i8,2(1x,f14.5)) go to 22 endif c c Now imode = 1 c else c c Set up the rotation matrices: c do is=1,nst(1) call setrot(ang1(is),ang2(is),ang3(is),anis1(is),anis2(is), + is,MAXROT,rotmat) end do c c Get ready to go: c call cova3(0.0,0.0,0.0,0.0,0.0,0.0,1,nst,MAXNST,c0,it,cc,aa, + 1,MAXROT,rotmat,cmax,maxcov) open(lout,file=outfl,status='UNKNOWN') c c Set the variogram data, direction by direction up to ndir directions: c i = 0 do id=1,ndir xx = -xoff(id) yy = -yoff(id) zz = -zoff(id) do il=1,nlag xx = xx + xoff(id) yy = yy + yoff(id) zz = zz + zoff(id) call cova3(0.0,0.0,0.0,xx,yy,zz,1,nst,MAXNST,c0,it, + cc,aa,1,MAXROT,rotmat,cmax,cov) i = i + 1 g(i) = dble(maxcov-cov) r(i) = dble(cov/maxcov) h(i) = sqrt(xx*xx+yy*yy+zz*zz) end do end do c c Now, loop over all the cutoffs: c if(icalc.eq.2) ncut=1 do i=1,ncut pcut = dble(zc(i)) if(icalc.eq.2) then pcut = dble(zcut) call coeffs(pcut,a,NTERMS) endif do l=1,ndir if(icalc.eq.1) then call dgauinv(pcut,xp,ierr) z = real(xp) write(lout,100) z,l else write(lout,101) zc(i),l endif 100 format('Model Indicator Variogram: at cutoff = ',f8.3, + ' Direction ',i2) 101 format('Model Normal Scores Variogram: for cutoff = ',f8.3, + ' Direction ',i2) do j=1,nlag ii = (l-1)*nlag + j if(icalc.eq.1) then rnew(ii) = rhoi(r(j),pcut,NTERMS) cnew(ii) = rnew(ii)*pcut*(1.0d0-pcut) gnew(ii) = (pcut*(1.0d0-pcut))-cnew(ii) gnewst(ii) = 1.0d0-rnew(ii) else rnew(ii) = rhoy(r(j),pcut,a,NTERMS) gnew(ii) = 1.0d0-rnew(ii) cnew(ii) = rnew(ii) gnewst(ii) = gnew(ii) endif if(icalc.eq.1.and.corr.eq..true.) then write(lout,102) j,h(ii),gnewst(ii),l, + rnew(ii),gnew(ii) else write(lout,102) j,h(ii),gnew(ii),l, + rnew(ii),gnewst(ii) endif 102 format(i3,1x,f8.3,1x,f12.3,1x,i6,4(1x,f12.5)) end do end do end do c c End of Main if statement c endif c c Finished: c close(lout) 9997 write(*,9998) VERSION 9998 format(/' BIGAUS2 Version: ',f5.3, ' Finished'/) stop 98 stop 'Error in parameter file somewhere' end subroutine coeffs(p,a,nterms) c--------------------------------------------------------------------------- c c Coefficients for Power Series Expansion c *************************************** c c Determines the coefficients for the power series approximation for c rhoi(ry,p) to be subsequently reverted (by subroutine rhoy(ri,p) c in order to obtain a power series for w=sqrt((1-ry)/(1+ry)). c c c INPUT PARAMETERS: c c p Cutoff/proportion c nterms Number of terms to include in series expansion c c c OUTPUT VARIABLES: c c a Coefficients for power series approximation: c the a(k) values are the coefficients in the c series for w as a function of the variable: c y = (kI*p*q+qq*qq-qq) c c c NOTES: c c The c(k) are the power series coefficients for rhoI in terms of w. c The a2k1 series is a polynomial approximation for atan(w), see c Abramowitz and Stegun (formula 4.4.49) c c c NO EXTERNAL REFERENCES c c Original: Marshal Grant 1996 c--------------------------------------------------------------------------- implicit double precision (a-h,o-z) parameter (ntmax=101) double precision ki,ki0 dimension c(ntmax),a2k1(17) dimension cmatrx(ntmax,ntmax) dimension a(nterms) data a2k1 / 1.0000000000d0, 0.0000000000d0,-0.3333314528d0, + 0.0000000000d0, 0.1999355085d0, 0.0000000000d0, + -0.1420889944d0, 0.0000000000d0, 0.1065626393d0, + 0.0000000000d0,-0.0752896400d0, 0.0000000000d0, + 0.0429096138d0, 0.0000000000d0,-0.0161657367d0, + 0.0000000000d0, 0.0028662257d0/ if (nterms.gt.ntmax) then write(6,*) 'In subroutine COEFFS, number of terms exceeds ntmax.' write(6,*) 'Increase ntmax = NTERMS' stop endif pi=4.0d0*datan(1.0d0) q=1.0d0-p qq=min(p,q) call dgauinv(p,xp,ierr) y0=xp bsq=0.5d0*y0*y0 b=sqrt(bsq) kk=1 term=1.0d0 correc=term csign=1.0d0 c(kk)=csign*(1.0d0-correc*dexp(-bsq))/(kk)-a2k1(kk) c(kk)=c(kk)/pi do 20 k=2,(nterms+1)/2 kk=2*k-1 c(kk-1)=0.0d0 term=term*bsq/(k-1) correc=correc+term csign=-csign c(kk)=csign*(1.0d0-correc*dexp(-bsq))/(kk) if (kk.le.17) c(kk)=c(kk)-a2k1(kk) c(kk)=c(kk)/pi 20 continue c c now do the reversion of series to get rhoY from rhoI c initialize cmatrx first c do 100 i=1,nterms do 110 j=1,nterms cmatrx(i,j)=0.0d0 110 continue 100 continue c c now load coefficients into first column c do 120 i=1,nterms cmatrx(i,1)=c(i) 120 continue do 200 jcol=2,nterms do 210 i1=1,nterms-1 do 220 i2=1,nterms-i1 irow=i1+i2 cmatrx(irow,jcol)= cmatrx(irow,jcol) + * c(i1)*cmatrx(i2,jcol-1) 220 continue 210 continue 200 continue c c now do the reversion c jcol=1 a(1)=1.0d0/cmatrx(1,1) do 300 jcol=2,nterms sum=0.0d0 do 310 i=1,jcol-1 sum=sum+a(i)*cmatrx(jcol,i) 310 continue a(jcol)=-sum/cmatrx(jcol,jcol) 300 continue return end double precision function rhoi(ry,p,nterms) c--------------------------------------------------------------------------- c c Indicator Correlation Coefficient c ********************************* c c Computes the indicator correlation coefficient for a bivariate standard c normal distribution. It does not actually develop a formal series c approximation but determines higher order terms by a recursion formula. c Although the number of terms (power of w) is given by nterms, the c routine actually computes only the (nterms+1)/2 nonzero terms. c c c INPUT PARAMETERS: c c ry Normal scores correlation coefficient c p Cutoff/proportion c nterms Number of terms to include in series expansion c c c c NO EXTERNAL REFERENCES c c Original: Marshal Grant 1996 c--------------------------------------------------------------------------- implicit double precision (a-h,o-z) double precision ki,ki0 nonzer=(nterms+1)/2 pi=4.0d0*datan(1.0d0) q=1.0d0-p qq=min(p,q) call dgauinv(p,xp,ierr) y0=xp bsq=0.5d0*y0*y0 b=sqrt(bsq) ki0=qq factor=dexp(-bsq) factor=(1.0d0/pi)*dexp(-bsq) wsq=(1.0d0-ry)/(1.0d0+ry) w=sqrt(wsq) fklast=atan(w) f0=fklast k=0 coeff=-factor ki=ki0+coeff*fklast do 10 k=1,nonzer fk=(w**(2*k-1))/(2*k-1)-fklast coeff=-coeff*bsq/k ki=ki+coeff*fk fklast=fk 10 continue rhoi=(ki-qq*qq)/(p*q) return end double precision function rhoy(ri,p,a,nterms) c-------------------------------------------------------------------------- c c Normal Scores Correlation Coefficient c ************************************* c c Computes the normal scores correlation coefficient, for a biavariate c Gaussian distribution, from the input cutoff/proportion and indicator c correlation coefficient. c c INPUT PARAMETERS: c c ri Indicator correlation coefficient c p Cutoff/proportion c a Vector of coefficients of series expansion c nterms Number of terms to include in series expansion c c c NO EXTERNAL REFERENCES c c Original: Marshal Grant 1996 c-------------------------------------------------------------------------- implicit double precision (a-h,o-z) dimension a(nterms) q=1.0d0-p qq=min(p,q) y=ri*p*q+qq*qq-qq w=0.0d0 do 10 i=nterms,1,-1 w=y*(a(i)+w) 10 continue wsq=w*w rhoy=(1.0d0-wsq)/(1.0d0+wsq) return end subroutine setrot(ang1,ang2,ang3,anis1,anis2,ind,MAXROT,rotmat) c----------------------------------------------------------------------- c c Sets up an Anisotropic Rotation Matrix c ************************************** c c Sets up the matrix to transform cartesian coordinates to coordinates c accounting for angles and anisotropy (see manual for a detailed c definition): c c c INPUT PARAMETERS: c c ang1 Azimuth angle for principal direction c ang2 Dip angle for principal direction c ang3 Third rotation angle c anis1 First anisotropy ratio c anis2 Second anisotropy ratio c ind matrix indicator to initialize c MAXROT maximum number of rotation matrices dimensioned c rotmat rotation matrices c c c NO EXTERNAL REFERENCES c c c----------------------------------------------------------------------- parameter(DEG2RAD=3.141592654/180.0,EPSLON=1.e-20) real*8 rotmat(MAXROT,3,3),afac1,afac2,sina,sinb,sint, + cosa,cosb,cost c c Converts the input angles to three angles which make more c mathematical sense: c c alpha angle between the major axis of anisotropy and the c E-W axis. Note: Counter clockwise is positive. c beta angle between major axis and the horizontal plane. c (The dip of the ellipsoid measured positive down) c theta Angle of rotation of minor axis about the major axis c of the ellipsoid. c if(ang1.ge.0.0.and.ang1.lt.270.0) then alpha = (90.0 - ang1) * DEG2RAD else alpha = (450.0 - ang1) * DEG2RAD endif beta = -1.0 * ang2 * DEG2RAD theta = ang3 * DEG2RAD c c Get the required sines and cosines: c sina = dble(sin(alpha)) sinb = dble(sin(beta)) sint = dble(sin(theta)) cosa = dble(cos(alpha)) cosb = dble(cos(beta)) cost = dble(cos(theta)) c c Construct the rotation matrix in the required memory: c afac1 = 1.0 / dble(max(anis1,EPSLON)) afac2 = 1.0 / dble(max(anis2,EPSLON)) rotmat(ind,1,1) = (cosb * cosa) rotmat(ind,1,2) = (cosb * sina) rotmat(ind,1,3) = (-sinb) rotmat(ind,2,1) = afac1*(-cost*sina + sint*sinb*cosa) rotmat(ind,2,2) = afac1*(cost*cosa + sint*sinb*sina) rotmat(ind,2,3) = afac1*( sint * cosb) rotmat(ind,3,1) = afac2*(sint*sina + cost*sinb*cosa) rotmat(ind,3,2) = afac2*(-sint*cosa + cost*sinb*sina) rotmat(ind,3,3) = afac2*(cost * cosb) c c Return to calling program: c return end real*8 function sqdist(x1,y1,z1,x2,y2,z2,ind,MAXROT,rotmat) c----------------------------------------------------------------------- c c Squared Anisotropic Distance Calculation Given Matrix Indicator c *************************************************************** c c This routine calculates the anisotropic distance between two points c given the coordinates of each point and a definition of the c anisotropy. c c c INPUT VARIABLES: c c x1,y1,z1 Coordinates of first point c x2,y2,z2 Coordinates of second point c ind The rotation matrix to use c MAXROT The maximum number of rotation matrices dimensioned c rotmat The rotation matrices c c c c OUTPUT VARIABLES: c c sqdist The squared distance accounting for the anisotropy c and the rotation of coordinates (if any). c c c NO EXTERNAL REFERENCES c c c----------------------------------------------------------------------- real*8 rotmat(MAXROT,3,3),cont,dx,dy,dz c c Compute component distance vectors and the squared distance: c dx = dble(x1 - x2) dy = dble(y1 - y2) dz = dble(z1 - z2) sqdist = 0.0 do i=1,3 cont = rotmat(ind,i,1) * dx + + rotmat(ind,i,2) * dy + + rotmat(ind,i,3) * dz sqdist = sqdist + cont * cont end do return end subroutine cova3(x1,y1,z1,x2,y2,z2,ivarg,nst,MAXNST,c0,it,cc,aa, + irot,MAXROT,rotmat,cmax,cova) c----------------------------------------------------------------------- c c Covariance Between Two Points c ***************************** c c This subroutine calculated the covariance associated with a variogram c model specified by a nugget effect and nested varigoram structures. c The anisotropy definition can be different for each nested structure. c c c c INPUT VARIABLES: c c x1,y1,z1 coordinates of first point c x2,y2,z2 coordinates of second point c nst(ivarg) number of nested structures (maximum of 4) c ivarg variogram number (set to 1 unless doing cokriging c or indicator kriging) c MAXNST size of variogram parameter arrays c c0(ivarg) isotropic nugget constant c it(i) type of each nested structure: c 1. spherical model of range a; c 2. exponential model of parameter a; c i.e. practical range is 3a c 3. gaussian model of parameter a; c i.e. practical range is a*sqrt(3) c 4. power model of power a (a must be gt. 0 and c lt. 2). if linear model, a=1,c=slope. c 5. hole effect model c cc(i) multiplicative factor of each nested structure. c (sill-c0) for spherical, exponential,and gaussian c slope for linear model. c aa(i) parameter "a" of each nested structure. c irot index of the rotation matrix for the first nested c structure (the second nested structure will use c irot+1, the third irot+2, and so on) c MAXROT size of rotation matrix arrays c rotmat rotation matrices c c c OUTPUT VARIABLES: c c cmax maximum covariance c cova covariance between (x1,y1,z1) and (x2,y2,z2) c c c c EXTERNAL REFERENCES: sqdist computes anisotropic squared distance c rotmat computes rotation matrix for distance c----------------------------------------------------------------------- parameter(PI=3.14159265,PMX=999.,EPSLON=1.e-10) integer nst(*),it(*) real c0(*),cc(*),aa(*) real*8 rotmat(MAXROT,3,3),hsqd,sqdist c c Calculate the maximum covariance value (used for zero distances and c for power model covariance): c istart = 1 + (ivarg-1)*MAXNST cmax = c0(ivarg) do is=1,nst(ivarg) ist = istart + is - 1 if(it(ist).eq.4) then cmax = cmax + PMX else cmax = cmax + cc(ist) endif end do c c Check for "zero" distance, return with cmax if so: c hsqd = sqdist(x1,y1,z1,x2,y2,z2,irot,MAXROT,rotmat) if(real(hsqd).lt.EPSLON) then cova = cmax return endif c c Loop over all the structures: c cova = 0.0 do is=1,nst(ivarg) ist = istart + is - 1 c c Compute the appropriate distance: c if(ist.ne.1) then ir = min((irot+is-1),MAXROT) hsqd=sqdist(x1,y1,z1,x2,y2,z2,ir,MAXROT,rotmat) end if h = real(dsqrt(hsqd)) c c Spherical Variogram Model? c if(it(ist).eq.1) then hr = h/aa(ist) if(hr.lt.1.) cova=cova+cc(ist)*(1.-hr*(1.5-.5*hr*hr)) c c Exponential Variogram Model? c else if(it(ist).eq.2) then cova = cova + cc(ist)*exp(-3.0*h/aa(ist)) c c Gaussian Variogram Model? c else if(it(ist).eq.3) then cova = cova + cc(ist)*exp(-(3.0*h/aa(ist)) + *(3.0*h/aa(ist))) c c Power Variogram Model? c else if(it(ist).eq.4) then cova = cova + cmax - cc(ist)*(h**aa(ist)) c c Hole Effect Model? c else if(it(ist).eq.5) then c d = 10.0 * aa(ist) c cova = cova + cc(ist)*exp(-3.0*h/d)*cos(h/aa(ist)*PI) cova = cova + cc(ist)*cos(h/aa(ist)*PI) endif end do c c Finished: c return end subroutine dgauinv(p,xp,ierr) c----------------------------------------------------------------------- c c Standard Normal Quantile c ************************ c c Computes the inverse of the standard normal cumulative distribution c function with a numerical approximation from : Statistical Computing, c by W.J. Kennedy, Jr. and James E. Gentle, 1980, p. 95. c c c INPUT/OUTPUT: c c p = double precision cumulative probability value: dble(psingle) c xp = G^-1 (p) in double precision c ierr = 1 - then error situation (p out of range), 0 - OK c c c----------------------------------------------------------------------- real*8 p0,p1,p2,p3,p4,q0,q1,q2,q3,q4,y,pp,lim,p,xp save p0,p1,p2,p3,p4,q0,q1,q2,q3,q4,lim c c Coefficients of approximation: c data lim/1.0e-20/ data p0/-0.322232431088/,p1/-1.0/,p2/-0.342242088547/, + p3/-0.0204231210245/,p4/-0.0000453642210148/ data q0/0.0993484626060/,q1/0.588581570495/,q2/0.531103462366/, + q3/0.103537752850/,q4/0.0038560700634/ c c Initialize: c ierr = 1 xp = 0.0d0 pp = p if(p.gt.0.5d0) pp = 1 - pp if(p.lt.lim) return ierr = 0 if(p.eq.0.5d0) return c c Approximate the function: c y = dsqrt(dlog(1.0/(pp*pp))) xp = y + ((((y*p4+p3)*y+p2)*y+p1)*y+p0) / + ((((y*q4+q3)*y+q2)*y+q1)*y+q0) if(p.lt.0.5d0) xp = -xp c c Return with G^-1(p): c return end subroutine chknam(str,len) c----------------------------------------------------------------------- c c Check for a Valid File Name c *************************** c c This subroutine takes the character string "str" of length "len" and c removes all leading blanks and blanks out all characters after the c first blank found in the string (leading blanks are removed first). c c c c----------------------------------------------------------------------- parameter (MAXLEN=132) character str(MAXLEN)*1 c c Remove leading blanks: c do i=1,len-1 if(str(i).ne.' ') then if(i.eq.1) go to 1 do j=1,len-i+1 k = j + i - 1 str(j) = str(k) end do do j=len,len-i+2,-1 str(j) = ' ' end do go to 1 end if end do 1 continue c c Find first blank and blank out the remaining characters: c do i=1,len-1 if(str(i).eq.' ') then do j=i+1,len str(j) = ' ' end do go to 2 end if end do 2 continue c c Return with modified file name: c return end subroutine makepar c----------------------------------------------------------------------- c c Write a Parameter File c ********************** c c c c----------------------------------------------------------------------- lun = 99 open(lun,file='bigaus2.par',status='UNKNOWN') write(lun,10) 10 format(' Parameters for BIGAUS2',/, + ' **********************',/,/, + 'START OF PARAMETERS:') write(lun,11) 11 format('1 ', + '-input mode') write(lun,12) 12 format('samplevar.out ', + '-file for input variogram') write(lun,13) 13 format('1 ', + '-threshold/proportion') write(lun,14) 14 format('1 ', + '-calculation mode') write(lun,15) 15 format('bigaus2.out ', + '-file for output of variograms') write(lun,16) 16 format('3 ', + '-number of thresholds') write(lun,17) 17 format('0.25 0.50 0.75 ', + '-threshold cdf values') write(lun,18) 18 format('1 20 ', + '-number of directions and lags') write(lun,19) 19 format('0.0 0.0 1.0 ', + '-azm(1), dip(1), lag(1)') write(lun,20) 20 format('2 0.2 ', + '-nst, nugget effect') write(lun,21) 21 format('1 0.4 0.0 0.0 0.0 ', + '-it,cc,ang1,ang2,ang3') write(lun,22) 22 format(' 10.0 5.0 10.0 ', + '-a_hmax, a_hmin, a_vert') write(lun,23) 23 format('1 0.4 0.0 0.0 0.0 ', + '-it,cc,ang1,ang2,ang3') write(lun,24) 24 format(' 10.0 5.0 10.0 ', + '-a_hmax, a_hmin, a_vert') close(lun) return end