c     gfortran pex_m4i_8ch_OP_heat.f m4iset_8ch.c /usr/lib64/libspcm_linux.so -o pex_m4i_OP_heat

      use iso_c_binding

      implicit real*8 (a-h,o-z)

      integer*2 ipex(10000000)
      integer*2 ipex1(10000000)
      integer*4 ndp
      
      real*8 amp(8),phi(8), phaseSub(8)
      complex*16 uvec(8), uvec0(8),uvecwc(8),chelp
      
      character*64 argu
      
      pi=4.d0*datan(1.d0)
      
      ipex=0
      ipex1=0
      
c      fsamp=625.d6
c      fsamp=500.d6
      fsamp=400.d6

c      fbas=123.246d6
      fbas=297.234d6

      ndp=512000
c      ndp=2*51200
      ndp8=ndp/8


      
      
      call getarg(1,argu)
      read(argu,*)mode


      call getarg(2,argu)
      read(argu,*)att


      amp0=att*32767.0d0

c     Calibration is off
c       phi(1) = phi(1) + 0.0
c       phi(2) = phi(2) + ((3.5)  * (pi / 180.0))
c       phi(3) = phi(3) + ((-7.0045896053314207)  * (pi / 180.0))
c       phi(4) = phi(4) + ((1.0142267853021623)   * (pi / 180.0))
c       phi(5) = phi(5) + ((24.745329284667967)   * (pi / 180.0))
c       phi(6) = phi(6) + ((16.720711231231689)   * (pi / 180.0))
c       phi(7) = phi(7) + ((10.251145744323731)   * (pi / 180.0))
c       phi(8) = phi(8) + ((8.3969257831573483)   * (pi / 180.0))

      if(mode.eq.1) then
      open(3,file='wcvec.dat',status='unknown')
      do i=1,8
      read(3,*)uvec_r,uvec_i
      uvecwc(i)=dcmplx(uvec_r,uvec_i)
      phaseSub(i)= uvec_i
      uvec0(i)=dsqrt(0.125d0)*cdexp((0.d0,-1.d0)*dfloat(i-1)*pi/4.d0)
      enddo
      
      
      close(3)




      do i=1,8
      chelp=(0.d0,0.d0)
      do ii=1,8
      chelp=chelp+conjg(uvecwc(ii))*uvec0(ii)
      enddo
c     WC/IF      
      uvec(i)=uvec0(i)-chelp*uvecwc(i)
c      uvec(i)=uvecwc(i)
      enddo
      
c     Calculate total Orthogonal Projection power 
      help1_OP=0.d0
      do i=1,8
      help1_OP=help1_OP+cdabs(uvec(i))**2
      enddo
      do i=1,8
      uvec(i)=uvec(i)/dsqrt(help1_OP)
      enddo 

      do i=1,8
c     CP
c      amp(i)=cdabs(uvec0(i))
c      phi(i)=datan2(dimag(uvec0(i)),dreal(uvec0(i)))
c     WC
c      amp(i)=dsqrt(0.125d0)*1.0
c      amp(i)=cdabs(uvecwc(i))
c      phi(i)=datan2(dimag(uvecwc(i)),dreal(uvecwc(i)))

C     OP 
      amp(i)=cdabs(uvec(i))
      phi(i)=datan2(dimag(uvec(i)),dreal(uvec(i)))
      
c      print *,amp(i)
c      print *,phi(i)


      enddo

      pow0=0.d0
      powwc=0.d0
      powmin=0.d0
      powWC_ONES = 0.d0
      do i=1,8
      pow0=pow0+cdabs(uvec0(i))**2
      powwc=powwc+cdabs(uvecwc(i))**2
      powmin=powmin+cdabs(uvec(i))**2
   
c      uvecwc = dcmplx(amp,phaseSub)
c      powWC_ONES=powWC_ONES+cdabs(uvecwc(i))**2
      powWC_ONES=powWC_ONES+amp(i)**2
      enddo
c      print*,uvecwc
      print*,pow0,powwc,powmin,powWC_ONES

      
      do i=1,ndp
      t=dfloat(i-1)/fsamp
      help1=dsin(2.d0*pi*fbas*t+phi(1))
      help2=dsin(2.d0*pi*fbas*t+phi(2))
      help3=dsin(2.d0*pi*fbas*t+phi(3))
      help4=dsin(2.d0*pi*fbas*t+phi(4))
      help5=dsin(2.d0*pi*fbas*t+phi(5))
      help6=dsin(2.d0*pi*fbas*t+phi(6))
      help7=dsin(2.d0*pi*fbas*t+phi(7))
      help8=dsin(2.d0*pi*fbas*t+phi(8))
      
      
c      ipex(4*i-3)=int(amp(1)*amp0*help1*(1.0/1.0))
c      ipex(4*i-2)=int(amp(2)*amp0*help2*(1.0/1.0))
c      ipex(4*i-1)=int(amp(3)*amp0*help3*(1.026722925))
c      ipex(4*i-0)=int(amp(4)*amp0*help4*(0.972435079))
c      ipex1(4*i-3)=int(amp(5)*amp0*help5*(0.970286109))
c      ipex1(4*i-2)=int(amp(6)*amp0*help6*(0.946317249))
c      ipex1(4*i-1)=int(amp(7)*amp0*help7*(1.006437117))
c      ipex1(4*i-0)=int(amp(8)*amp0*help8*(0.961280215))

      ipex(4*i-3)=int(amp(1)*amp0*help1)
      ipex(4*i-2)=int(amp(2)*amp0*help2)
      ipex(4*i-1)=int(amp(3)*amp0*help3)
      ipex(4*i-0)=int(amp(4)*amp0*help4)
      ipex1(4*i-3)=int(amp(5)*amp0*help5)
      ipex1(4*i-2)=int(amp(6)*amp0*help6)
      ipex1(4*i-1)=int(amp(7)*amp0*help7)
      ipex1(4*i-0)=int(amp(8)*amp0*help8)
      enddo

      else

      do i=1,ndp
      t=dfloat(i-1)/fsamp
c      help=dsin(2.d0*pi*fbas*t)
      help1=dsin(2.d0*pi*fbas*t+phi(1))
      help2=dsin(2.d0*pi*fbas*t+phi(2))
      help3=dsin(2.d0*pi*fbas*t+phi(3))
      help4=dsin(2.d0*pi*fbas*t+phi(4))
      help5=dsin(2.d0*pi*fbas*t+phi(5))
      help6=dsin(2.d0*pi*fbas*t+phi(6))
      help7=dsin(2.d0*pi*fbas*t+phi(7))
      help8=dsin(2.d0*pi*fbas*t+phi(8))
c      if(i/ndp8.eq.0)ipex(4*i-3)=int(amp0*help1*(1.0))
c      if(i/ndp8.eq.1)ipex(4*i-2)=int(amp0*help2*(0.9776674937))
c      if(i/ndp8.eq.2)ipex(4*i-1)=int(amp0*help3*(1.0287206266))
c      if(i/ndp8.eq.3)ipex(4*i-0)=int(amp0*help4*(0.93213295074))
c      if(i/ndp8.eq.4)ipex1(4*i-3)=int(amp0*help5*(0.533))
c      if(i/ndp8.eq.5)ipex1(4*i-2)=int(amp0*help6*(0.95037037))
c      if(i/ndp8.eq.6)ipex1(4*i-1)=int(amp0*help7*(0.95039))
c      if(i/ndp8.eq.7)ipex1(4*i-0)=int(amp0*help8*(0.935))

      if(i/ndp8.eq.0)ipex(4*i-3)=int(amp0*help1)
      if(i/ndp8.eq.1)ipex(4*i-2)=int(amp0*help2)
      if(i/ndp8.eq.2)ipex(4*i-1)=int(amp0*help3)
      if(i/ndp8.eq.3)ipex(4*i-0)=int(amp0*help4)
      if(i/ndp8.eq.4)ipex1(4*i-3)=int(amp0*help5)
      if(i/ndp8.eq.5)ipex1(4*i-2)=int(amp0*help6)
      if(i/ndp8.eq.6)ipex1(4*i-1)=int(amp0*help7)
      if(i/ndp8.eq.7)ipex1(4*i-0)=int(amp0*help8)
      enddo
      
      endif


      do i=1,1
      print *,i,'----------------------------'
      aread=m4iset(ipex,ipex1,ndp)
      print *,i,'----------------------------'
      enddo

      end

