Next: REFERENCES Up: Code Appendix Previous: Makefile

Variogram code

```Program variogram

use sep
use bin1
use globecalc

implicit none
real                              :: az1,az2,ov,dv,a,b,c,d,pi
real,dimension (:)  , pointer     :: lat,lon,shaz,temp,anglesum,ang,delta,ones,fold
integer                           :: idat,idat2,nv,n1,k,stat

call sep_init ()

call from_history("n1",n1)

allocate (lat(n1),lon(n1),shaz(n1))

!write(0,*) 'lat is ', lat
!write(0,*) 'lon is ', lon
!write(0,*) 'az is  ', shaz

call from_param( "nv", nv)
call from_param( "ov", ov)
call from_param( "dv", dv)

call to_history("n1",nv)
call to_history("n2", 1)
call to_history("o1",ov)
call to_history("d1",dv)

call to_history("label2","orthogonality")
call to_history("label1","kilometers")

call sep_close()
pi=3.141593
!------------quadrant tests for the subtraction ------
!  k=1
!  if (k/=0) then
!     a=asin(sin( 0.17453295-pi )) !10 in radians
!     b=asin(sin( 3.3161259 -pi )) !190 in radians
!     a=asin(sin( 10*(pi/180) ))                 ! = 0.17453295
!     b=asin(sin(-10*(pi/180) ))                 ! =-0.17453295
!     c=asin(sin(170*(pi/180) ))                 ! = 0.174532562
!     d=asin(sin(190*(pi/180) ))                 ! =-0.174533218
!     write(0,*) '10 -90=',a*(180/pi),' 190 -90=',b*(180/pi),'   170=',c,'   190=',d
!
!  else   !---to use, must also take out end if at end of program--
!-----------------------------------

lat= lat * (pi/180)             !---Convert degrees lat/long to radians --
lon= lon * (pi/180)
shaz=shaz* (pi/180)

write(0,*) 'Doing the variance bit'
allocate(delta ( n1*(n1-1)/2) )
allocate(ang   ( n1*(n1-1)/2) )
delta=0.  ; ang=0.
!     write(0,*) 'allocated'
!     write(0,*) 'size of ang is ',size(ang),'delta is ',size(delta)
k=0
do idat =1,n1-1
!        write(0,*) 'Finished ',k,' of',size(ang),' loops.'
if( shaz(idat) <=180 ) then                 !null values in data file =999
do idat2=idat+1,n1
k=k+1
if( shaz(idat2) <=180 ) then
call calc( lat(idat),lon(idat),lat(idat2),lon(idat2),delta( (idat-1)*n1+(idat2-1) ),az1,az2)
ang(k) = sqrt( abs(asin(sin( shaz(idat)-shaz(idat2) ))))
end if
end do
end if
end do
!---------------------------------------------
!---------------putting data into distance bins ---------------
delta= delta * 6371                             !  convert to km's: rad. of earth=6371 km's
allocate( anglesum(nv) )
anglesum=0                                      !  initialize to zero
!     write(0,*)'size of anglesum is ', size(anglesum)
call bin1_init( nv, ov, dv, delta )             !  bins should go to 0.5*circ. of earth
stat = bin1_lop( .true., .true.,anglesum, ang)  !  anglesum is output
!--------------to normalize anlgesum,  fold=Btrans (ones) -----------
allocate( ones(n1) )
allocate( fold(nv) )
ones=1.; fold=0.
call bin1_init( nv, ov, dv, delta )
stat = bin1_lop( .true., .false., fold, ones)
!  write(0,*) 'anglesum at 1 is ', anglesum
where(fold/=0)                                 !  mind not dividing by zero
anglesum=anglesum/fold                      !  divide by number of data point in the bin
elsewhere
anglesum=0
end where
!  write(0,*) 'fold at 1 is ', fold
!  write(0,*) 'anglesum is now', anglesum

call sep_write( anglesum)
!  end if   !---for when you don't really want to run the program
end program variogram
```

Next: REFERENCES Up: Code Appendix Previous: Makefile
Stanford Exploration Project
3/24/2001