program lineofsightyinyang_fs

!==========================================================================
! Program to generate integrated lines of sight through a 3-D netcdf data
! file. Reads data in from file, calculates the line of sight using the
! ionization timescale of ejecta, then creates a 2-D netcdf image containing
! the output.

! Note also: lineofsight assumes that 1) the data to be read is is in
! spherical geometry, 2) the y and z coordinates correspond to polar and
! azimuthal angles respectively, 3) the netCDF dataset has a fixed-width
! grid, and 4) the data is typical yin-yang output.

!==========================================================================

IMPLICIT NONE

logical :: onYINgrid

character(len=4) :: outtimestep, timestep, tmp1
character(len=50) :: prefix, geomfile, casefile, varname
character(LEN=80) :: label1a, label1b, label1c, label1d, label1e, label1f, label1g, label1h
character(LEN=80) :: label2, label3, label4
character(len=16) :: coord1, coord2, coord3, coord4
character(len=19) :: ghostfilename, varfilename
character(len=22) :: infilename

integer :: parts
integer :: n, i, ks, ktmp
integer :: ncid, ncstat
integer :: ndim, nvar, nv, rank, nfile, nstart, nstop
integer :: imax, jmax, kmax, nmax, ipart, mfile
integer, dimension(:,:), allocatable :: ongrid

real, dimension(:,:,:), allocatable :: yin, yang
real(8), dimension(:), allocatable :: zxc, zyc, zzc

! variables for line of sight integration
integer :: dbg_int ! DEBUG LINE
integer :: ig, jg, kg, j, k, is, missed
integer :: ib, ia, jb, ja, kb, ka

real :: pi, temp1, temp2, alpha, delta, pcutoff, g1, g2, rmax, time, expand
real :: xmin, xmax, ymin, ymax, delx, dely, delz, rz, rad, phi, the, phi_yin, the_yin
real :: radmin, radmax, themin, themax, phimin, phimax, xy, yy, zy, tx, ty, tz, xn, yn, zn
real :: c1, c2, c3, c4, c5, c6, r1, r2, r3, r4, t1, t2
real, dimension(:), allocatable :: x, y
real, dimension(:,:), allocatable :: emission

! variables for emission data writeout
character(len=12) :: vname
character(len=25) :: outfilename

pi = 2.0 * asin(1.0)
!=========================================================================
write(6,*) 'Create a series of line-of-sight images'
write(6,fmt="('Input PREFIX of EnSight data file: (without leading E or F)')")
read (5,*) prefix

write(6,fmt="('Input data file number (from zero): ')")
read (5,*) nfile
write(timestep,"(i4)") nfile + 1000

!infilename = trim(prefix) // '.rmax'
!open(17,file=infilename)
!read(17,*) mfile, time, rmax
!do while (mfile < nfile+1000)
!  read(17,*) mfile, time, rmax
!enddo
!print *, 'using rmax = ', rmax, ' from file ',mfile
!close(17)

!=========================================================================

! set size (physical and pixel) of images
ig   = 1500  ! resolution of final image
jg   = 1500  ! resolution of final image

! set rotation angle and change in angle between frames
alpha = 0.0
delta = 0.0*pi / 160.

!==========================================================================
! OPEN GHOST FILE AND DETERMINE ONGRID STATUS
!==========================================================================

! OPEN DATA FILE with ghost zones and geometry
ghostfilename = 'Laurenghost.bin'
print *, 'opening ', ghostfilename
open(unit=15,file=ghostfilename,form='unformatted')
  read(15) imax, jmax, kmax, ks

  ! allocate coordinate arrays
  ALLOCATE( zxc(imax) )
  ALLOCATE( zyc(jmax) )
  ALLOCATE( zzc(kmax) )
  ! allocate data array for variable
  ALLOCATE(ongrid(jmax,kmax) )
 
  read(15) zxc, zyc, zzc
  read(15) ongrid
close(15)

print *, imax, jmax, kmax

!expand = rmax / zxc(imax)
!do i = 1, imax
! zxc(i) = zxc(i) * expand
!enddo
rmax = zxc(imax)

xmax = rmax * 1.1 
ymax = xmax
xmin = -xmax
ymin = -ymax

ALLOCATE( x(ig) )
ALLOCATE( y(jg) )
ALLOCATE( emission(ig,jg) )

delx = (xmax-xmin)/real(ig)
dely = (ymax-ymin)/real(jg)
do i = 1, ig 
  x(i) = xmin + real(i)*delx
enddo
do j = 1, jg 
  y(j) = ymin + real(j)*dely
enddo
!==========================================================================
  ! allocate arrays
  ALLOCATE (  yin(imax,jmax,kmax) )
  ALLOCATE ( yang(imax,jmax,kmax) )

  varfilename = 'E' // trim(prefix) // timestep // '.rho'
  open(unit=15,file=varfilename,form='unformatted')
   read(15) varname
   read(15) label2
   read(15) ipart
   read(15) label3
   read(15) yin
  close(15)

  varfilename = 'F' // trim(prefix) // timestep // '.rho'
  open(unit=15,file=varfilename,form='unformatted')
   read(15) varname
   read(15) label2
   read(15) ipart
   read(15) label3
   read(15) yang
  close(15)



!========================================================================
! PERFORM LINE OF SIGHT INTEGRATION
!========================================================================

print *, 'integrate along line of sight'


! calculate zone widths to be used in linear interpolations later
delx  = zxc(2) - zxc(1)
dely  = zyc(2) - zyc(1)
delz  = zzc(2) - zzc(1)

! define limits of simulation grid
radmax = zxc(imax-1)
radmin = zxc(2)
themax = zyc(jmax-1)
themin = zyc(2)
phimax = zzc(kmax-1)
phimin = zzc(2)

emission = 0.0

! Start main integration loop, taking care that when an index is 0 angles
!   and trig functions don't explode.

  do i = 1, ig
    ! With loops nested i/k/j, lineofsight looks down the polar axis in 
    ! its integration. If loops are nested i/j/k, integration is along an
    ! equatorial axis.  To switch between the two, change nesting of loops
    ! and change indices of emission in both locations where it is used.

!$OMP PARALLEL DO PRIVATE(j,k,tx,ty,tz,xn,yn,zn,rz,rad,phi,the,c1,c2,c3,c4,c5,c6,t1,t2,ib,ia,jb,ja,kb,ka,onYINgrid)

    do j = 1, jg
    ! emission should be using i,j indices

      do k = 1, ig   ! loop along line of sight

!#########################################################################
        ! apply any coordinate transformations (e.g. rotations) here if
        ! transforms apply to cartesian coordinates
        tx = x(i)
        ty = y(j)
        tz = x(k)

        xn = +tx * cos(alpha) + tz *sin(alpha)
        yn = +ty 
        zn = -tx * sin(alpha) + tz * cos(alpha)
!#########################################################################

        ! calculate location of cell <x(i),x(j),x(k)> in spherical coords
        rz  = sqrt(xn**2 +         zn**2)
        rad = sqrt(xn**2 + yn**2 + zn**2)
        if (xn /= 0) then
          phi = atan2(zn,xn) ! phi returned in (-pi, pi)
        else
          if (zn <  0.) phi = -pi*0.5 ! x(i)=0, x(k)<0
          if (zn >= 0.) phi =  pi*0.5 ! x(i)=0, x(k)>0
        endif
        if (rz /= 0.) then
          the = atan(yn/rz) ! theta returned in (-pi/2,pi/2)
          the = 0.5*pi - the  ! map (pi/2,-pi/2) to (0,pi)
        else
          if (yn < 0.) the = pi
          if (yn >= 0.) the = 0.
        endif

        ! is integration on one of the simulation grids?
        if( (rad > radmin) .and. (rad < radmax) ) then

          ib = 1 + (rad-zxc(1))/delx ! zone below <x(i),x(j),x(k)> in radius
          ia = 1 + ib                ! zone above in radius
          c1 = (zxc(ia) - rad)/delx
          c2 = (rad - zxc(ib))/delx

          ! is it on the yin (sunny) or yang grid?  requires two conditions, inside the angle limits, and inside ongrid
          onYINgrid = .false.
          if ((the>themin).and.(the<themax).and.(phi>phimin).and.(phi<phimax)) then 
            jb = 1 + (the-zyc(1))/dely ! zone below cell in theta
            kb = 1 + (phi-zzc(1))/delz ! zone below in phi
            ja = 1 + jb                ! zone above in theta
            ka = 1 + kb                ! zone above in phi

            ! if zones bordering <x(i),x(j),x(k)> are not ghost zones
            if( ongrid(jb,kb)==0 .and. ongrid(jb,ka)==0 .and. &
                ongrid(ja,kb)==0 .and. ongrid(ja,ka)==0 )        onYINgrid = .true.
          endif
              

          if (onYINgrid) then

             ! calculate coefficients for linear interpolation
             c3 = (zyc(ja) - the)/dely
             c4 = (the - zyc(jb))/dely
             c5 = (zzc(ka) - phi)/delz
             c6 = (phi - zzc(kb))/delz

             ! perform the linear interpolation to estimate emission variable
             t1 = c3*(yin(ib,jb,kb)*c1 + yin(ia,jb,kb)*c2) + c4*(yin(ib,ja,kb)*c1 + yin(ia,ja,kb)*c2)
             t2 = c3*(yin(ib,jb,ka)*c1 + yin(ia,jb,ka)*c2) + c4*(yin(ib,ja,ka)*c1 + yin(ia,ja,ka)*c2)
             emission(i,j) = emission(i,j) + (c5*t1 + c6*t2)**2

          else ! it must be on the yang grid

!#########################################################################
           ! transform theta and phi to the yang grid
           !   (apply any coordinate transformations (e.g. rotations) here as
           !   well)
           tx = -x(i)   ! x -> -x
           ty = x(k)    ! y -> z
           tz = y(j)    ! z -> y

           xn = tx * cos(alpha) - ty * sin(alpha)
           yn = tx * sin(alpha) + ty * cos(alpha)
           zn = tz
!#########################################################################

           rz  = sqrt(xn**2+zn**2)
           if (xn /= 0.) then
             phi = atan2(zn,xn) ! phi returned in (-pi, pi)
           else
             if (zn <  0.) phi = -pi*0.5 ! x(i)=0, x(k)<0
             if (zn >= 0.) phi =  pi*0.5 ! x(i)=0, x(k)>0
           endif
           if (rz /= 0.) then
             the = atan(yn/rz) ! theta returned in (-pi/2,pi/2)
             the = 0.5*pi - the  ! map (pi/2,-pi/2) to (0,pi)
           else
             if (yn <  0.) the = pi
             if (yn >= 0.) the = 0.
           endif

           if ((the>themin).and.(the<themax).and.(phi>phimin).and.(phi<phimax)) then

           jb = 1 + (the-zyc(1))/dely ! zone below cell in theta
           kb = 1 + (phi-zzc(1))/delz ! zone below in phi
           ja = 1 + jb                ! zone above in theta
           ka = 1 + kb                ! zone above in phi

           ! if zones bordering <x(i),x(j),x(k)> are not ghost zones
           if( ongrid(jb,kb)==0 .and. ongrid(jb,ka)==0 .and. &
               ongrid(ja,ka)==0 .and. ongrid(ja,ka)==0 ) then

             ! calculate coefficients for linear interpolation
             c3 = (zyc(ja) - the)/dely
             c4 = (the - zyc(jb))/dely
             c5 = (zzc(ka) - phi)/delz
             c6 = (phi - zzc(kb))/delz

             ! perform the linear interpolation to estimate emission variable
             t1 = c3*(yang(ib,jb,kb)*c1 + yang(ia,jb,kb)*c2) + c4*(yang(ib,ja,kb)*c1 + yang(ia,ja,kb)*c2)
             t2 = c3*(yang(ib,jb,ka)*c1 + yang(ia,jb,ka)*c2) + c4*(yang(ib,ja,ka)*c1 + yang(ia,ja,ka)*c2)
             emission(i,j) = emission(i,j) + (c5*t1 + c6*t2)**2
 
           else
             print *, 'not ongrid for either yin or yang'
           endif ! ongrid? check

           else
             print *, 'fell through the baseball seam'
           endif

          endif ! on the yin grid

        endif ! in netcdf file? check

      enddo

    enddo
!$OMP END PARALLEL DO 

  enddo 

alpha = alpha + delta

!==========================================================================
! WRITE EMISSION TO NEW EnSight GOLD FILE
!==========================================================================
ipart = 1
label3 = 'block'
label2 = 'part'
label4 = 'Emission'

outfilename = 'T' // trim(prefix) // timestep // '.los'
open(unit=15,file=outfilename,form='unformatted')
 write(15) label4
 write(15) label2
 write(15) ipart
 write(15) label3
 write(15) emission
close(15)

!==========================================================================
! Write an Ensight geometry file

! set up labels for EnSight files
label1a = 'Fortran Binary'
label1b = 'cartesian image file from los code'
label1c = 'junk line 2'
label1d = 'node id off'
label1e = 'element id off'
label1f = 'part '
label1g = 'cartesian coordinates'
label1h = 'block rectilinear'

kg = 1
xn = 0.0
! open up the ensight geometry file
geomfile = 'T' // trim(prefix) // '.geo'
open(unit=15,file=geomfile,form='unformatted')
  write(15) label1a
  write(15) label1b
  write(15) label1c
  write(15) label1d
  write(15) label1e
  write(15) label1f
  write(15) ipart
  write(15) label1g
  write(15) label1h
  write(15) ig, ig, kg
  write(15) x
  write(15) x
  write(15) xn
close(15)

! Write an Ensight case file

casefile = 'T' // trim(prefix) // '.case'
open(unit=14,file=casefile)
  write(14,*) 'FORMAT'
  write(14,*) 'type: ensight gold'
  write(14,*) 'GEOMETRY'
  write(14,*) 'model: ',  trim(geomfile) 
  write(14,*) 'VARIABLE'
  write(14,*) 'scalar per node:    emission   ', trim(outfilename)
close(14)

!#############################################################


stop

end program lineofsightyinyang_fs