      program main
c This is modified from distfn.f.
      parameter(npx=513)
      implicit real * 8 (a-h,o-z)
      common/refpara/a1,a2,a3,a4,a5,a6,a7,ijk
      common / comalt / alt(npx,npx),xp(npx),yp(npx)
      dimension blt(5,npx,npx)
      common /pda/x1min,xxmin,xxmax,y1min,yymin,yymax,nx,ny
c     read(5,*) xmin,xmax,ymin,ymax,nx,ny
c     xxmin=140
c     xxmax=145
c     yymin=35
c     yymax=40
c     xxmin=128
c     xxmax=149
c     yymin=30
c     yymax=47
c     xxmin=0
c     xxmax=360
c     yymin=-90
c     yymax=90
c     nx=100
c     ny=100
      do 100 ijk=1,1
c     do 100 ijk=1,5
c     do 100 ijk=3,3
      call datfij(tx,ty,alt,ijk,nd,npex)
      fmin=100
      fmax=0
      do i=1,nx
      do j=1,ny
      blt(ijk,i,j)=alt(i,j)
      if(fmin.gt.alt(i,j)) fmin=alt(i,j)
      if(fmax.lt.alt(i,j)) fmax=alt(i,j)
      end do
      end do
      write(6,*) 'fmin',fmin,fmax
c     if(ijk.eq.1) write(11,*) nx,ny
c     write(11,2)((alt(i,j),i=1,nx),j=1,ny)
    2 format(4e20.10)
  100 continue
      open(10,file='interpolated.pixel')
      do i=1,nx
      do j=1,ny
      write(10,1) xp(i)+x1min,yp(j)+y1min,(blt(k,i,j),k=1,1)
    1 format(2f9.3,5e13.5)
      end do
      end do
      close(10)
      stop
      end
      subroutine datfij(tx,ty,fij,ijk,nd,npex)
c  this is function on the delaunay(delone) triangulation
c  of 2-dimensional space. (12/march/89)
c  estimated parameter data is inoue.data
c  parameter data is inotest.data
c     program fdelon
      implicit real * 8 (a-h,o-z)
      parameter(mdl=244444,npa=323456,npx=513)
      common/refpara/a1,a2,a3,a4,a5,a6,a7,ijk0
      dimension ldt(mdl,3),adt(mdl),idt(mdl)
      common /pda/x1min,xxmin,xxmax,y1min,yymin,yymax,nx,ny
      common / comalt / alt(npx,npx),xp(npx),yp(npx)
      dimension xx(npa),yy(npa),zxy(npa),zxy1(npa)
      dimension xi(npx),yj(npx),fij(npx,npx)
c
c     nx=100
c     ny=100
      if(ijk.eq.1)
     &call redata(np,npex,nd,ldt,adt,idt,
     &                               xx,yy,zxy1,tx,ty)
c     call redata(np,npex,nd,ldt,adt,idt,xx,yy,zxy1,tx,ty)
c
      zmin=100
      zmax=0
      do 10 i=1,npex
      zxy(i)=zxy1(i+(ijk-1)*npex)
      if(zmin.gt.zxy(i)) zmin=zxy(i)
      if(zmax.lt.zxy(i)) zmax=zxy(i)
   10 continue
c     write(6,*) 'i',zxy(1),zxy(npex),npex
      write(6,*) ijk,zmin,zmax
      call flattc(tx,ty,nd,ldt,xx,yy,zxy,fij,xp,yp)
      return
      end
c*********************************************************************
      subroutine redata(np,npex,nd,ldt,adt,idt,
     &                               xx,yy,zxy,tx,ty)
c  inputs:
c         np: number of points within a rectangular
c       npex: number of points including the minimum # of mirror image.
c         nd: number of the delaunay triangulars.
c         tx: horizontal length of the rectangular.
c         ty: vertical lengh of the rectangular.
c       dens: density=np/(tx*ty)
c        ipt: identification of a triangular: i=ipt in the following.
c   ldt(i,.): identification of the verteces of the i-th triangular.
c   ndt(i,.): id of the neighboring triangular of the i-th one.
c     adt(i): area of the i-th triangular
c     mdt(j); # of the triangulars including the vertex j.
c     idt(j): id of the triangulars using the vertex j.
c
      implicit real * 8 (a-h,o-z)
      parameter(mdl=244444,npa=323456)
      common /pda/x1min,xxmin,xxmax,y1min,yymin,yymax,nx,ny
      common/refpara/a1,a2,a3,a4,a5,a6,a7,ijk
      dimension ldt(mdl,3),ndt(mdl,3),adt(mdl),mdt(mdl),idt(mdl)
      dimension xx(npa),yy(npa),zz(npa,5),zxy(npa),sd(npa),mp(npa)
      character*50 file1,file2
c reading configulation of the delaunay tassilation
      open(8,file='interpolated.conf')
      read(8,1) file1
      read(8,1) file2
      read(8,*) x1min,y1min,xxmin,xxmax,yymin,yymax,nx,ny
      write(6,2) 'x,y,nx,ny',x1min,y1min,xxmin,xxmax,yymin,yymax,nx,ny
    2 format(a,6f10.2,2i5)
c     if(xxmin.lt.0) xxmin=360+xxmin
c     if(xxmax.lt.0) xxmax=360+xxmax
c     yymin=yymin+90
c     yymax=yymax+90
c     x1min=120
c     y1min=20
      xxmin=xxmin-x1min
      xxmax=xxmax-x1min
      yymin=yymin-y1min
      yymax=yymax-y1min
    1 format(a)
      close(8)
c     open(4,file='alljpM50n.config')
      open(4,file=file1)
      read(4,*) np,npex,nd,tx,ty
      read(4,*) (idum,xx(ipt),yy(ipt),d1,d2,d3,d4,d5,ipt=1,npex)
      do i=1,nd
      read(4,*) idum,(ldt(i,j),j=1,3),adt(i)
      enddo
      close(4)
c
c     open(9,file='delopara.test')
      open(9,file=file2)
      read(9,*)  w1,abic,mx
      write(6,*)  'w1',w1,abic,mx
c     read(9,*)  w1,w2,w3,w4,w5,abic,mx
c     read(9,*) a1,a2,a3,a4,a5,a6,a7
c     read(9,*) a1,a2,a3,a4,a5,a6,a7
c     read(9,*) a1,a2,a3,a4,a5,a6,a7
c     read(9,*) a1,a2,a3,a4,a5,a6,a7
c     read(9,*) a1,a2,a3,a4,a5,a6,a7
      read(9,*) (zxy(i), i=1,mx)
      close(9)
c
      write(6,4) tx,ty,npex,nd
      write(6,*)  w1,abic,mx
c     write(6,*) a1,a2,a3,a4,a5,a6,a7
    4 format(1h ,'tx,ty =',2f8.1,'  npex =',i10,' nd',i10)
c     do 10 i=1,npex
c       if(mod(i,100).eq.0) write(6,5) i,xx(i),yy(i),zxy(i)
c  10 continue
    5 format(1h ,i7,3f10.3)
      return
      end
c*********************************************************************
      subroutine flattc(tx,ty,nd,ldt,xx,yy,zxy,fij,xp,yp)
c-----------------------------------------------------------------------
c      calculating values on the given lattice points
c      for drawing contour maps or bird's-eye-view.
c-----------------------------------------------------------------------
c  inputs:
c       npex: number of points including the minimum # of mirror image.
c         nd: number of the delaunay triangulars.
c         tx: horizontal length of the rectangular.
c         ty: vertical lengh of the rectangular.
c       dens: density=np/(tx*ty)
c        ipt: identification of a triangular: i=ipt in the following.
c   ldt(i,.): identification of the verteces of the i-th triangular.
c   ndt(i,.): id of the neighboring triangular of the i-th one.
c
c   output
c     fij(i,j): values of the delaunay surface on the lattices (xi,yj)
c
      implicit real * 8 (a-h,o-z)
      parameter(mdl=244444,npa=323456,npx=513)
      common/refpara/a1,a2,a3,a4,a5,a6,a7,ijk
      common /pda/x1min,xxmin,xxmax,y1min,yymin,yymax,nx,ny
      dimension ldt(mdl,3)
      dimension xx(npa),yy(npa),zxy(npa)
      dimension fij(npx,npx),xp(npx),yp(npx)
      write(6,*) xxmin,xxmax
      write(6,*) yymin,yymax
c     stop
      do i=1,nx
c     xii=tx*(i-0.5d0)/nx
      xii=xxmin+(xxmax-xxmin)*(i-0.5d0)/nx
c     xii=xxmin+(xxmax-xxmin)*(i-1)/(nx-1)
      xp(i)=xii
      do j=1,ny
c     yjj=ty*(j-0.5d0)/ny
      yjj=yymin+(yymax-yymin)*(j-0.5d0)/ny
c     yjj=yymin+(yymax-yymin)*(j-1)/(ny-1)
      yp(j)=yjj
c     write(6,*) i,j,xii,yjj
c     if(ijk.eq.3)write(6,*) i,j,xii,yjj
      call interp(nd,xx,yy,zxy,ldt,xii,yjj,fijij)
      fij(i,j)=fijij
c     write(6,*) i,j,xii,yjj,fij(i,j),a1*exp(fij(i,j))
   55 format('fij',3f8.2)
      end do
      end do
      do 111 i=1,nx
      do 111 j=1,ny
      if(ijk.eq.1) fij(i,j)=fij(i,j)
c     if(ijk.eq.1) fij(i,j)=a1*exp(fij(i,j))
c     if(ijk.eq.1) fij(i,j)=a1*exp(fij(i,j))
c     if(ijk.eq.2) fij(i,j)=a2*exp(fij(i,j))
c     if(ijk.eq.3) fij(i,j)=a4*exp(fij(i,j))
c     if(ijk.eq.4) fij(i,j)=a5*exp(fij(i,j))
c     if(ijk.eq.5) fij(i,j)=a7*exp(fij(i,j))
  111 continue
      return
      end
      subroutine interp(nd,xx,yy,zxy,ldt,xi,yj,fij)
c  inputs:
c         nd: number of the delaunay triangulars.
c   xx(i)
c   yy(i)
c   zxy(i)
c   ldt(i,.): identification of the verteces of the i-th triangular.
c   xi 
c   yj
c  output
c     fij: values of the delaunay surface on the lattices (xi,yj)
c
      implicit real * 8 (a-h,o-z)
      parameter(mdl=244444,npa=323456)
      dimension ldt(mdl,3)
      dimension xx(npa),yy(npa),zxy(npa)
      dimension xd(3),yd(3),zd(3)
      det(x1,x2,x3,y1,y2,y3)=x1*y2+x2*y3+x3*y1-x1*y3-x2*y1-x3*y2
      data eps/-0.1d-7/
      do ipt=1,nd
        do k=1,3
          xd(k)=xx(ldt(ipt,k))
          yd(k)=yy(ldt(ipt,k))
          zd(k)=zxy(ldt(ipt,k))
        end do
        denom=det(xd(1),xd(2),xd(3),yd(1),yd(2),yd(3))
      if(denom.eq.0.0) write(6,*) 'denominator=0.0 for ipt =', ipt
      if(denom.eq.0.0) go to 20
        aa1=det(xi,xd(2),xd(3),yj,yd(2),yd(3))/denom
        aa2=det(xd(1),xi,xd(3),yd(1),yj,yd(3))/denom
        aa3=det(xd(1),xd(2),xi,yd(1),yd(2),yj)/denom
c           if(a1.lt.0. .or. a2.lt.0. .or. a3.lt.0.) go to 20
      if(aa1.lt.eps.or. aa2.lt.eps.or. aa3.lt.eps) go to 20
        fij=aa1*zd(1)+aa2*zd(2)+aa3*zd(3)
c      write(6,2) ipt,aa1,aa2,aa3,fij,zd(1),zd(2),zd(3)
c      write(6,*) xi,yj
    2 format('ipt',i5,7f10.2,e20.10)
c      write(6,*) (ldt(ipt,k),k=1,3)
c      write(6,*) (xx(ldt(ipt,k)),k=1,3)
c      write(6,*) (yy(ldt(ipt,k)),k=1,3)
c      write(6,*) (zxy(ldt(ipt,k)),k=1,3)
        return
   20 continue
      end do
      write(6,*) 'error'
      stop 119
      end
