      program pgen
*     ------- ------

      include 'Online.cmm'

      character*15 tem
      character*12 ras
      integer system
	common /esrwateg/specor,ynom(200),xnom(200)

        common /pawc/ hmem(500000)
        call hlimit(500000)
	

         open(34,file='special_cor.txt')
           read(34,*)specor
         close(34)

	write(*,*)'Generate Deformation Plots '
        if (specor.ne.0.) then
	write(*,*)'All RASAS Y are corrected with ',specor,
     +                     ' um!!!!!!'
	endif
	


         open(34,file='timelast.txt')
           read(34,*)timet
         close(34)

	call surfview(1,timet)
	call surfview(2,timet)
	call surfview(3,timet)
	call surfview(4,timet)

      close(35)
      lun  = 35
      CALL HROPEN(lun,'oops','pgen.uit','N',1024,ISTAT)
      if (istat.ne.0) then
        CALL HROPEN(lun,'oops','pgen.uit','U',1024,ISTAT)
      endif
      call hrout(0,icycle,'T')
      call hrend('oops')
      close(lun)

      end

	subroutine checknom(igroup)
	integer ncheck,idcheck(200)
        ncheck=0
        do j=1,4
	  k=(j-1)*4
	  if (j.eq.4) k=20
        do i=1,4
	  idcheck(i)=k+i ! RASAS
          ncheck=ncheck+1
        enddo
        enddo

	do i=37,45
	  idcheck(i-32)=i ! InPlane and XP rasniks
          ncheck=ncheck+1
	enddo

	do i=1,ncheck
           write(*,*)'id= ',id',' dy ',bla blas


	
	

	  
	subroutine surfview(igroup,timef)
        include 'Online.cmm'
	character*80 fnaam,fnaam1,fnaam2,fnaam3
        character*500 line
        character*240 dline,sline
        character*3 iline
        common /esrwateg/specor,ynom(200),xnom(200)
        real xras(0:3)
        real yras(0:3)
        real zras(0:3)
        real rotZras(0:3)
        integer icol(51)
	integer iap(4)
	integer iam(4)
	integer icm(4)
	integer icp(4)
        integer irasas(4)
        data iap/ 1,2,3,4/
        data iam/ 5,6,7,8/
        data icm/ 9,10,11,12/
        data icp/ 21,22,23,24/
* 33 34 35 36
* 37 38 39 40
* 41 42 43 44
* 13 14 15 16
* 17 18 19 20
* 25 26 27 28
* 29 30 31 32
	real ylast(200),xlast(200),slast(200),dtime(200),
     +       rlast(200),dtfralast(200)
	time=timef
	do i=1,200
	  dtime(i)=1000.
	  ynom(i)=0.
	  xnom(i)=0.
	  ylast(i)=999.
	  xlast(i)=0.
	  rlast(i)=0.
	  slast(i)=0.
	enddo

	open(39,file='nominal.dat')
 14	continue
	read(39,*,ERR=15,END=15)ig,id,y,x
           if (ig.eq.igroup .or. ig.eq.0) then
	     ynom(id)=y
	     xnom(id)=x
           endif
        goto 14
 15	continue
	close(39)
             xnom(37)=20.
             xnom(38)=21.
             xnom(39)=-79.
             xnom(40)=20.

	if (igroup.eq.1) then
	     fnaam2='inarrow1.kumac'
	     fnaam1='surfcol1.kumac'
	     fnaam3='inpcol1.kumac'
	     fnaam='surfarrow1.kumac'
             do i=1,4
		irasas(i)=iap(i)
             enddo
             ynom(37)=-70.
             ynom(38)=-71.
             ynom(39)=-79.
             ynom(40)=-52.
	elseif (igroup.eq.2) then
	     fnaam2='inarrow2.kumac'
	     fnaam1='surfcol2.kumac'
	     fnaam3='inpcol2.kumac'
	     fnaam='surfarrow2.kumac'
             do i=1,4
                irasas(i)=iam(i)
             enddo
             ynom(37)=60.
             ynom(38)=74.
             ynom(39)=76.
             ynom(40)=57.
        elseif (igroup.eq.3) then
	     fnaam2='inarrow3.kumac'
	     fnaam1='surfcol3.kumac'
	     fnaam3='inpcol3.kumac'
	     fnaam='surfarrow3.kumac'
             do i=1,4
                irasas(i)=icm(i)
             enddo
             ynom(37)=60.
             ynom(38)=74.
             ynom(39)=76.
             ynom(40)=67.
        elseif (igroup.eq.4) then
	     fnaam2='inarrow4.kumac'
	     fnaam1='surfcol4.kumac'
	     fnaam3='inpcol4.kumac'
	     fnaam='surfarrow4.kumac'
             do i=1,4
                irasas(i)=icp(i)
             enddo
             ynom(37)=-50.
             ynom(38)=-76.
             ynom(39)=-74.
             ynom(40)=-72.
        else 
             print *,'no group selected in surfview'
	     stop
	endif

	print *,' Prepare plot for group ',igroup
	print *,' Time requested ',time

C--
C first fill rasnik vars from ntuple
C--
        ii=0
        close(35)
        nfound=0
        
        open(35,file='ntuple.ascii')
 10     continue
                read(35,'(a)',ERR=20,END=20)line
                write(iline,*)line(2:4)
                write(dline,*)line(10:249)
                read(iline,*,ERR=20,END=20)idch
                read(dline,*,ERR=20,END=20)
     +        timent,dt,
     +        dtgra,dtfra,
     +        xp,yp,scale,rotz,igch,ilay
               ii=ii+1
             if (abs(time-timent).lt.dtime(idch)) then ! take it
C 45=MILO
*                print *,'about to take ',timent,' for id ',idch
*                print *,'for ig igold ',igch,igroup_old
                if (igroup.eq.igch .or. 
     +              (idch.eq.45 .and. igroup_old.eq.igroup)) then
                  print *,'took time ',timent,' for id ',idch
		  xlast(idch)=xp
		  ylast(idch)=yp
		  slast(idch)=scale
		  rlast(idch)=rotz
                  dtfralast(idch)=dtfra
                  dtime(idch)=abs(time-timent)
                  nfound=nfound+1
                endif ! igch
            endif ! lowest time
	    igroup_old=igch
        goto 10
 20     continue
        close(35)

	ynom(1)=6.
	ynom(2)=6.
	ynom(3)=6.
	ynom(4)=6.
        ynom(5)=-6.
        ynom(6)=-6.
        ynom(7)=-6.
        ynom(8)=-6.
        ynom(9)=-6.
        ynom(10)=-6.
        ynom(11)=-6.
        ynom(12)=-6.
        ynom(21)=6.
        ynom(22)=6.
        ynom(23)=6.
        ynom(24)=6.

	do ii=1,200
	    ynom(ii)=ynom(ii)/1000.
            if (xnom.ne.999) then
	       xnom(ii)=xnom(ii)/1000.
           endif
	enddo




	nzbin=20
	nxbin=20
	yvalmin=10000.
	yvalmax=-10000.
	call hbook1(igroup+60,'surface',nxbin,0.,5.,0.)
	call hbook1(igroup+70,'surface',nzbin,0.,2.,0.)
	call hbook2(igroup+50,'surface',nxbin,0.,5.,
     +                                  nzbin,0.,2.,0.)
	call hbook2(igroup+40,'surface',nxbin,0.,5.,
     +                                  nzbin,0.,2.,0.)
        call hbook2(igroup+30,'inplane',nxbin,0.,5.,
     +                                  nzbin,0.,2.,0.)

        do iz=1,nzbin
           do ix=1,nxbin
             x=(ix-1)*5.0/((nxbin-1)*1.)
             z=(iz-1)*2.0/((nzbin-1)*1.)
             dxx=x/5.0
             dzz=z/2.0
             c=xlast(irasas(4)) ! SW
             a=xlast(irasas(3)) ! SE
             d=xlast(irasas(1)) ! NW
             b=xlast(irasas(2)) ! NE
             
             call fpint(xval,dxx,dzz,a,b,c,d)
	     c=ylast(irasas(4)) ! SW
	     a=ylast(irasas(3)) ! SE
	     d=ylast(irasas(1)) ! NW
	     b=ylast(irasas(2)) ! NE
             call fpint(yval,dxx,dzz,a,b,c,d)
             rdiv=1.
	     if (ylast(44).ne.0 .and. ylast(43).ne.0) rdiv=2.
             ysaghv=(ylast(44)+1.5*ylast(43))/rdiv/2. ! image shift
             rdiv=1.
	     if (ylast(42).ne.0 .and. ylast(41).ne.0) rdiv=2.
             ysagro=(ylast(42)+1.5*ylast(41))/rdiv/2. ! image shift
C MILO 45
	     ysagmi=0.5*ylast(45)*(1030./800.)**2
             ysagmi=ysagmi*(1.-(1.-z)**2)
             ysagro=ysagro*(1.-(1.-z)**2)
             ysaghv=ysaghv*(1.-(1.-z)**2)
             if (igroup.eq.2 .or. igroup.eq.3) then
                ysagmi=-ysagmi
                ysagro=-ysagro
		ysaghv=-ysaghv
	     endif
             ysagz=ysagro
	     ysagn=ysaghv
             if (igroup.eq.2 .or. igroup.eq.4) then
		ysagz=ysaghv
		ysagn=ysagro
	     endif
             if (dxx.gt.0.5) then
                dxxx=dxx*2.-1.
                yval=yval+ysagn*(dxxx)**2
                yval=yval+ysagmi*(2.*(0.5-abs(dxx-0.5)))**2
             else
                dxxx=dxx*2.
                yval=yval+ysagz*(1.-dxxx)**2
                yval=yval+ysagmi*(2.*dxx)**2
             endif
             yval=yval*1000.
             call hix(60+igroup,ix,xf)
             call hix(70+igroup,iz,zf)
             xf=xf+0.01
             zf=zf+0.01
             yval=yval-specor ! special correction
	     call hfill(50+igroup,xf,zf,yval)
	     if (yval.gt.yvalmax) yvalmax=yval
	     if (yval.lt.yvalmin) yvalmin=yval
	     call hfill(40+igroup,xf,zf,xval*1000.)
	enddo
	enddo


        call prepal(fnaam1,50+igroup,6.,yvalmin,yvalmax)
*            1 ------- 2
*            |         |
*            3 ------- 4   function values
         
C prepare arrows to represent rasasX in PAW
        x1=0.
        y1=0.27

	x2=0.63
        y2=0.45

	x3=0.37
	y3=0.0015

	x4=1.
	y4=0.15

	
	open(35,file=fnaam)
        dxne=x2-x4
        dyne=y2-y4
        rn=sqrt(dxne**2+dyne**2)
        b=xlast(irasas(2)) ! NE
        write(35,*)'color 1 0 0 1 '
        if (b.gt.0.012) write(35,*)'color 1 1 0 0 '
        if (b.lt.0.006) write(35,*)'color 1 0 1 0 '
        b=b*rn/0.100
        x=x4-b*dxne/rn
        y=y4-b*dyne/rn
	write(35,*)'arrow ',x4,x,y4,y
C frascatie
        xf=(x4+x2)/2.
        yf=(y4+y2)/2.
        b=dtfralast(irasas(2))*23.*2.2/1000.
        write(35,*)'color 1 0 0 1 '
        if (b.gt.0.012) write(35,*)'color 1 1 0 0 '
        if (b.lt.0.006) write(35,*)'color 1 0 1 0 '
        b=b*rn/0.100
        x=xf-b*dxne/rn
        y=yf-b*dyne/rn
        write(35,*)'arrow ',xf,x,yf,y
C


        b=xlast(irasas(1)) ! NW
        write(35,*)'color 1 0 0 1 '
        if (b.gt.0.012) write(35,*)'color 1 1 0 0 '
        if (b.lt.0.006) write(35,*)'color 1 0 1 0 '
        b=b*rn/0.100
        x=x2-b*dxne/rn
        y=y2-b*dyne/rn
	write(35,*)'arrow ',x2,x,y2,y

        dxse=x1-x3
        dyse=y1-y3
        rn=sqrt(dxse**2+dyse**2)
        b=xlast(irasas(3)) ! SE
        write(35,*)'color 1 0 0 1 '
        if (b.gt.0.012) write(35,*)'color 1 1 0 0 '
        if (b.lt.0.006) write(35,*)'color 1 0 1 0 '
        b=b*rn/0.100
        x=x3-b*dxse/rn
        y=y3-b*dyse/rn
        write(35,*)'arrow ',x3,x,y3,y

C frascatie
        b=dtfralast(irasas(3))*23.*2.2/1000.
        write(35,*)'color 1 0 0 1 '
        if (b.gt.0.012) write(35,*)'color 1 1 0 0 '
        if (b.lt.0.006) write(35,*)'color 1 0 1 0 '
        b=b*rn/0.100
        xf=(x3+x1)/2.
        yf=(y3+y1)/2.
        x=xf-b*dxse/rn
        y=yf-b*dyse/rn
        write(35,*)'arrow ',xf,x,yf,y



        b=xlast(irasas(4)) ! SW
        write(35,*)'color 1 0 0 1 '
        if (b.gt.0.012) write(35,*)'color 1 1 0 0 '
        if (b.lt.0.006) write(35,*)'color 1 0 1 0 '
        b=b*rn/0.100
        x=x1-b*dxse/rn
        y=y1-b*dyse/rn
        write(35,*)'arrow ',x1,x,y1,y

	write(35,*)'color 1 0 0 0 '
	close(35)

C now the inplanestuff
        do i=0,3
	  xras(i)=xlast(37+i)
	  yras(i)=ylast(37+i)
          if (slast(37+i).ne.0.) then
            zras(i)=2500./slast(37+i)-2500.
          endif
          rotZras(i)=rlast(37+i)/1000. ! rotzrad in radians
        enddo
        icyc=1
        call get_deform_one(igroup,icyc,xras,yras,zras,rotZras)

          yvalmax=-100000.
          yvalmin=100000.

        do iz=1,nzbin
           do ix=1,nxbin
             x=(ix-1)*5.0/((nxbin-1)*1.)
             z=(iz-1)*2.0/((nzbin-1)*1.)
             xx=x-5.0/2.
             zz=z-2.0/2.
             yy=0.
             call shift(igroup,xx,yy,zz,xuse,yuse,zuse)
             call hix(60+igroup,ix,xf)
             call hix(70+igroup,iz,zf)
             xf=xf+0.01
             zf=zf+0.01
             yval=yuse*1000.
             call hfill(30+igroup,xf,zf,yval)
             if (yval.gt.yvalmax) yvalmax=yval
             if (yval.lt.yvalmin) yvalmin=yval
          enddo
          enddo

           call prepal(fnaam3,30+igroup,50.,yvalmin,yvalmax)

C prepare arrow for paw
	dbin=50./1000. ! in um  unless in prepal
	dtwo=2.*dbin/1000. ! in um  unless in prepal
        x1=0.
        y1=0.27
        x2=0.63
        y2=0.45
        x3=0.37
        y3=0.0015
        x4=1.
        y4=0.15
        open(35,file=fnaam2)
        dxne=x2-x4
        dyne=y2-y4
        rn=sqrt(dxne**2+dyne**2)
        b=chasgz(igroup,1)
        write(35,*)'color 1 0 0 1 '
        if (b.gt.dtwo) write(35,*)'color 1 1 0 0 '
        if (b.lt.dbin) write(35,*)'color 1 0 1 0 '
        b=b*rn/0.100
        xc=(x3+x4)/2.
        yc=(y3+y4)/2.
        x=xc-b*dxne/rn
        y=yc-b*dyne/rn
        write(35,*)'arrow ',xc,x,yc,y
        xc=(x1+x2)/2.
        yc=(y1+y2)/2.
        x=xc-b*dxne/rn
        y=yc-b*dyne/rn
        write(35,*)'arrow ',xc,x,yc,y

        b=temccd(igroup,1)*2.*0.025
        write(35,*)'color 1 0 0 1 '
        if (b.gt.dtwo) write(35,*)'color 1 1 0 0 '
        if (b.lt.dbin) write(35,*)'color 1 0 1 0 '
        b=b*rn/0.100
        x=x1-b*dxse/rn
        y=y1-b*dyse/rn
        write(35,*)'arrow ',x1,x,y1,y
        x=x3+b*dxse/rn
        y=y3+b*dyse/rn
        write(35,*)'arrow ',x3,x,y3,y



        b=temmsk(igroup,1)*2.*0.025
        write(35,*)'color 1 0 0 1 '
        if (b.gt.dtwo) write(35,*)'color 1 1 0 0 '
        if (b.lt.dbin) write(35,*)'color 1 0 1 0 '
        b=b*rn/0.100
        x=x4-b*dxse/rn
        y=y4-b*dyse/rn
        write(35,*)'arrow ',x4,x,y4,y
        x=x2+b*dxse/rn
        y=y2+b*dyse/rn
        write(35,*)'arrow ',x2,x,y2,y


	close(35)


	end







*-- Author :    Marcel Vreeswijk   03/03/94
      SUBROUTINE FPINT(fval,dxx,dyy,a,b,c,d)
******************************************************************
*                                                                *
* FPINT (version 1.00)                                           *
*                                                                *
*                                                                *
* IN:
* interpolate linear on a four point grid                        *
* dyy= dy/ybinwid
* dxx= dx/xbinwid
*            c ------- d
*            |         |
*            a ------- b   function values                       *
* OUT:   fval corresponding to dxx dyy                           *
*                                                                *
*                                                                *
*  written by: Marcel Vreeswijk, 03/03/94 10.53.52               *

      real Fval,dxx,dyy,a,b,c,d
      Fval  = a+ dxx*(b-a)+ dyy*(
     +     (c+ dxx * (d -c) )-
     +     (a+ dxx * (b -a)) )
      end




	include 'prepal.f'
	include 'deform_shift.f'
