      program main
      implicit real*8(a-h,o-z)
      parameter(npa=12345,nda=12345,mdl=33333,nr=030)
      parameter(nm=12345)
      common /coef/ aa(nda),bb(nda),cc(nda)
      common /domain/ amulog(nm),aklog(nm),alflog(nm),plog(nm),
     &                qlog(nm),ainteg,bgrmax,mx,my,nmain
      common /triang/ adt(mdl),ldt(mdl,3),ndt(mdl,3),nn,npext,ndelo
      common /basepa/ amu0,ak0,c,alf0,p0,d,q0
      common /orgepi/ alon(npa),alat(npa),amag(npa)
      common /delink/ nzu(npa),iu(nr,npa)
      dimension xd(3),yd(3),zm(3),zd(3),za(3),zp(3),zq(3),zb(3)
      dimension xx(npa),yy(npa),bv(npa),ldtb(mdl,3)
      dimension xxb(npa),yyb(npa)
      dimension xi(50000),xj(50000),bvxy(50000)
      det(x1,x2,x3,y1,y2,y3)=x1*y2+x2*y3+x3*y1-x1*y3-x2*y1-x3*y2
      eps=1.d-15
      open(9,file='hist-etas-mk-forecast.conf')
      read(9,*) tx,ty,tz
      read(9,*) xxmin,yymin
      close(9)
      nxy=0
      open(3,file='./hist-etas-mk.upda')
      read(3,*) w1,w2,abic, npa2
      read(3,*) amu0,ak0,c,alf,p,d,q
      read(3,*) (amulog(i),i=1,npa2/2),(aklog(i),i=1,npa2/2)
      close(3)
      open(1,file='./delone2.out')
      read(1,*) nn,npext,ndelo,tx,ty
      do 21 i=1,npext
      read(1,*) idum,alon(i),alat(i),amag(i),ot,dd,ee,ff
      if(i.gt.nn) go to 21
      aa(i)=ee/dd/sqrt(1-ff**2)
      bb(i)=ff   /sqrt(1-ff**2)
      cc(i)=dd/ee/sqrt(1-ff**2)
   21 continue
      do 509 i=1,ndelo
      read(1,*) idummy,(ldt(i,j),j=1,3),adt(i)
  509 continue
      do 520 i = 1, npext
        read(1,*) idummy, nzu(i), (iu(j,i), j=1,nzu(i))
  520 continue
      close(1)
      open(2,file='etas2aniso.out3')
      open(4,file='work.para')
      n1=0
      read(2,*)
  250 continue
      read(2,*,end=260)  idum,along,alatt,amg,xt,ax2,bx2,cx2
      n1=n1+1
      x=along-xxmin
      y=alatt-yymin
      do 220 ipt=1,ndelo
        xmin=1.d69
        ymin=1.d69
        xmax=-1.d69
        ymax=-1.d69
        do 230 k=1,3
          xd(k)=alon(ldt(ipt,k))
          yd(k)=alat(ldt(ipt,k))
          zd(k)=aklog(ldt(ipt,k))
          zm(k)=amulog(ldt(ipt,k))
          zb(k)=bv(ldt(ipt,k))
          if(xd(k).lt.xmin) xmin=xd(k)
          if(xd(k).gt.xmax) xmax=xd(k)
          if(yd(k).lt.ymin) ymin=yd(k)
          if(yd(k).gt.ymax) ymax=yd(k)
  230     continue
           if(x.lt.xmin) go to 220
           if(x.gt.xmax) go to 220
           if(y.lt.ymin) go to 220
           if(y.gt.ymax) go to 220
           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 220
           aa1=det(x    ,xd(2),xd(3),y    ,yd(2),yd(3))/denom
           aa2=det(xd(1),x    ,xd(3),yd(1),y    ,yd(3))/denom
           aa3=det(xd(1),xd(2),x    ,yd(1),yd(2),y    )/denom
           if(aa1.lt.eps.or. aa2.lt.eps.or. aa3.lt.eps) go to 220
           aklog1=aa1*zd(1)+aa2*zd(2)+aa3*zd(3)
           amulog1=aa1*zm(1)+aa2*zm(2)+aa3*zm(3)
           go to 240
  220    continue
  240 continue
      xt=xt+tz
      bv1=log(0.9d0)
      write(6,8) n1,along,alatt,amg,xt,ax2,bx2,cx2,
     &           amulog1,aklog1
    8 format(i5,f9.4,f8.4,f4.1,f12.5,3f9.4,2x,2e15.6)
      write(4,9) amulog1,aklog1,
     &           along,alatt,amg,xt,ax2,bx2,cx2
    9 format(2e15.6,2x,f9.4,f8.4,f4.1,f12.5,3f9.4)
      go to 250
  260 close(2)
      close(4)
    7 format(2f7.2,d20.12)
      stop
      end
