From d5d219b9c6e2267e6e35a890d2d3a5e00f833e9f Mon Sep 17 00:00:00 2001
From: Berk Silemek <berk.silemek@gmail.com>
Date: Mon, 14 Nov 2022 15:36:57 +0000
Subject: [PATCH] Fortran code to produce a single shot worst-case
 RF-excitation vector

---
 Software/Transmitter/pex_m4i_8ch_WC.f | 207 ++++++++++++++++++++++++++
 1 file changed, 207 insertions(+)
 create mode 100644 Software/Transmitter/pex_m4i_8ch_WC.f

diff --git a/Software/Transmitter/pex_m4i_8ch_WC.f b/Software/Transmitter/pex_m4i_8ch_WC.f
new file mode 100644
index 0000000..5885349
--- /dev/null
+++ b/Software/Transmitter/pex_m4i_8ch_WC.f
@@ -0,0 +1,207 @@
+c     gfortran pex_m4i_8ch_WC.f m4iset_8ch.c /usr/lib64/libspcm_linux.so -o pex_m4i_WC
+
+      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/OP      
+c      uvec(i)=uvec0(i)-chelp*uvecwc(i)
+      uvec(i)=uvecwc(i)
+      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
+      amp(i)=cdabs(uvecwc(i))
+      phi(i)=datan2(dimag(uvecwc(i)),dreal(uvecwc(i)))
+
+C     IF 
+c      amp(i)=cdabs(uvec(i))
+c      phi(i)=datan2(dimag(uvec(i)),dreal(uvec(i)))
+
+      print *,amp(i)
+      print *,phi(i)
+
+
+      enddo
+
+       
+c       phi(1) = phi(1) + 0.0
+c       phi(2) = phi(2) + ((3.25)  * (pi / 180.0))
+c       phi(3) = phi(3) + ((2.45)  * (pi / 180.0))
+c       phi(4) = phi(4) + ((0.0)   * (pi / 180.0))
+c       phi(5) = phi(5) + ((5.867324722)   * (pi / 180.0))
+c       phi(6) = phi(6) + ((5.167324722)   * (pi / 180.0))
+c       phi(7) = phi(7) + ((1.5)   * (pi / 180.0))
+c       phi(8) = phi(8) + ((1.5)   * (pi / 180.0))
+
+c       
+c      do i=1,8
+c      amp(i) = amp(i)*(dsqrt(0.84860632790893242d0))
+c      print *,amp(i)
+c      print *,phi(i)
+c      print*, uvecwc(i)
+c      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
+      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
+
-- 
GitLab