!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! Finite-Difference-Time -Domain Method
! Fortran Code
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

program pbg
use fgl
use fml
implicit none
!******************fortran data***************
complex(8) :: Ey(251,410)
complex(8) :: Hx(251,411),Hxi(251,2),Hxo(251,2)
complex(8) :: Hz(252,410),Hzi(2,410),Hzo(2,410)
complex(8) :: v
real(8) :: V1(251,410),V2(251,410)
real(8) :: r(251,410),E0(251,1),x(1,251), z(410)
real(8) :: a,cdt,dx,dz,Lx,Lz,P,RU,q,q0,k0,x0
real(8) :: ratio1,ratio2,xr
integer(4) :: IU(250),JU(409),i(251),j(410),mz(408)
integer(4) :: MU,NU,x1,Loop,n1,n2,n3,p1,p2,lp
logical(4) :: first, anim
real(8) :: t1, t2
!***********************************************
real(8) :: n(size(x),size(z))
integer(4) :: z1(20),len
type(mfArray) :: MFX,MFY,MF1,MF2,h1,h2,zr

call CPU_TIME(t1)
v=MF_I;
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! Input Parameters

q=1.3; !Incident wavelength(um)
k0=2*MF_PI/q;
q0=1.3;
ratio1=0.35;
ratio2=0.2;
a=ratio1*q0;
RU=ratio2*a;

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!Calculation x-z Domain and Grid Size

dx=0.02;
dz=0.02;
p1=11;
p2=18;

Lx=p1*a;
Lz=p2*a;

x=mfColon(-Lx/2,dx,Lx/2);
z=mfColon(0.0,dz,Lz);

MU=size(x)
NU=size(z)



i=mfColon(1,MU)
IU=mfColon(2,MU)
j=mfColon(1,NU)
JU=mfColon(2,NU)

mz=mfColon(3,NU)
cdt=dz/5;
P=(cdt-dz)/(cdt+dz);

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! Refractive Index Distribution PBG
! in air
n=1.d0
r=0.d0

! in silicon
do n1=1,p1
  do n2=1,p2
    do n3=0,180,5
      zr=mfColon(-RU*sin(n3*MF_PI/180),dz,RU*sin(n3*MF_PI/180));
      len = mfSize(zr)
      xr=RU*cos(n3*MF_PI/180);
      x1=mfRound(mf((n1*a-a/2+xr)/dx));
      z1(1:len)=mfRound((n2*a-a/2+zr)/dz);
      if (n1==4 .and. n2<=10) then
      elseif (n1>=4 .and. n2==10) then
      else
        n(x1,z1(1:len)) = sqrt(11.4)
      end if
    end do !n3
  end do !n2
end do !n2

! Initial Condition

V1=n**2+cdt*r/2*120*MF_PI;
V2=n**2-cdt*r/2*120*MF_PI;
Hz=0;
Hx=0;
Ey=0;
x0=2*a;
E0=transpose(exp(-((x+x0)/a)**2));

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! Propagation Begin
Loop=anint(Lz/cdt)*2;
first=.true.
anim=.true.

call msMeshGrid(mfOut(MFX,MFY),mf(z(mz)),mf(x))
call msAssign(MF1,n(:,mz))

! call msRecordStart('pbg.avi')
do lp=1,Loop
  ! 1st run, percent to complete
  ratio1=dble(lp)/dble(Loop)
  write(*,*) ratio1
  Hxi(i,1:2)=Hx(i,1:2);
  Hxo(i,1:2)=Hx(i,NU:NU+1);
  Hzi(1:2,j)=Hz(1:2,j);
  Hzo(1:2,j)=Hz(MU:MU+1,j);

  Ey(i,2)=Ey(i,2)+E0(:,1)*exp(v*k0*lp*cdt);
  Hx(i,JU)=Hx(i,JU)+cdt/dz*(Ey(i,JU)-Ey(i,JU-1));

  Hz(IU,j)=Hz(IU,j)-cdt/dx*(Ey(IU,j)-Ey(IU-1,j));
  Ey=cdt*((Hx(i,j+1)-Hx(i,j))/dz-(Hz(i+1,j)-Hz(i,j))/dx)/V1+Ey*V2/V1;

  Hx(i,1)=Hxi(i,2)+P*(Hx(i,2)-Hxi(i,1));
  Hx(i,NU+1)=Hxo(i,1)+P*(Hx(i,NU)-Hxo(i,2));
  Hz(1,j)=Hzi(2,j)+P*(Hz(2,j)-Hzi(1,j));
  Hz(MU+1,j)=Hzo(1,j)+P*(Hz(MU,j)-Hzo(2,j));

! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! % Display

  if(first) then
    call msAssign(MF2,real(Ey(:,mz)))
    h1=mfSolidContour(MFX,MFY,MF1);
    call msGset(h1,'iso',mf((/0.0,2.0/)))
    call msHold('on')
    h2=mfFastPColor(MF2, mf((/ 0.04, 8.18, -Lx/2, Lx/2 /)))
    call msColormap('gray')
    call msColormapRange(-2,3)
    call msColorbar('vert')
    call msAxis('off')
    call msDrawMaterial(h2,'edge','visible','off')
    call msCamZoom(3)
    call msDrawNow()
    first=.false.
  else
    if (anim .and. mod(lp,10)==0) then
      call msAssign(MF2,real(Ey(:,mz)))
      call msGSet(h2,'cdata',MF2)
      call msDrawNow()
    end if
  end if

end do  !end of Loop
call msAssign(MF2,real(Ey(:,3:NU)))
call msGSet(h2,'cdata',MF2)
call msDrawNow()
! call msRecordEnd()

call CPU_TIME(t2)
write (*,*) 't=', t2-t1
call msViewPause()

end program
