Merge branch 'options-for-reading-in-photon-modes-should-have-the-same-features' into 'develop'

external photon mode same as photon modes in input file

See merge request !1489

src/poisson/photon_mode.F90

@@ -125,20 +125,28 @@ contains
           if (ncols == 5) then
             read(iunit, *) this%omega(ii), this%lambda(ii), &
               this%pol(ii,1), this%pol(ii,2), this%pol(ii,3)
-
-            ! Normalize polarization vector
-            this%pol = this%pol/sqrt(dot_product(this%pol(ii,:), this%pol(ii,:)))
-
-            ! Calculate polarization dipole
-            do ip = 1, mesh%np
-              this%pol_dipole(ip, ii) = dot_product(this%pol(ii, 1:this%dim), mesh%x(ip, 1:this%dim))
-            end do
+          else if (ncols == 7) then
+            this%has_q0_p0 = .true.
+            SAFE_ALLOCATE(this%pt_coord_q0(1:this%nmodes))
+            SAFE_ALLOCATE(this%pt_momen_p0(1:this%nmodes))
+            read(iunit, *) this%omega(ii), this%lambda(ii), &
+              this%pol(ii,1), this%pol(ii,2), this%pol(ii,3), &
+              this%pt_coord_q0(ii), this%pt_momen_p0(ii)
           else
             ! error if not 5 columns
             message(1) = 'Error: unexpected number of columns in file:'
             message(2) = filename
             call messages_fatal(2)
           end if
+
+          ! Normalize polarization vector
+          this%pol = this%pol/sqrt(dot_product(this%pol(ii,:), this%pol(ii,:)))
+
+          ! Calculate polarization dipole
+          do ip = 1, mesh%np
+            this%pol_dipole(ip, ii) = dot_product(this%pol(ii, 1:this%dim), mesh%x(ip, 1:this%dim))
+          end do
+
         end do
         call io_close(iunit)
       end if
@@ -150,10 +158,20 @@ contains
         SAFE_ALLOCATE(this%omega(1:this%nmodes))
         SAFE_ALLOCATE(this%lambda(1:this%nmodes))
         SAFE_ALLOCATE(this%pol(1:this%nmodes,1:3))
+        SAFE_ALLOCATE(this%pol_dipole(1:mesh%np, 1:this%nmodes))
+        if (ncols == 7) then
+          SAFE_ALLOCATE(this%pt_coord_q0(1:this%nmodes))
+          SAFE_ALLOCATE(this%pt_momen_p0(1:this%nmodes))
+        end if
       end if
       call MPI_Bcast(this%omega(1), this%nmodes, MPI_FLOAT, 0, mpi_world%comm, ierr)
       call MPI_Bcast(this%lambda(1), this%nmodes, MPI_FLOAT, 0, mpi_world%comm, ierr)
       call MPI_Bcast(this%pol(1,1), this%nmodes*3, MPI_FLOAT, 0, mpi_world%comm, ierr)
+      call MPI_Bcast(this%pol_dipole, mesh%np*this%nmodes, MPI_FLOAT, 0, mpi_world%comm, ierr)
+      if (ncols == 7) then
+        call MPI_Bcast(this%pt_coord_q0(1), this%nmodes, MPI_FLOAT, 0, mpi_world%comm, ierr)
+        call MPI_Bcast(this%pt_momen_p0(1), this%nmodes, MPI_FLOAT, 0, mpi_world%comm, ierr)
+      end if
 #endif
     else
       if (.not. parse_is_defined(namespace, 'PhotonModes')) then
@@ -194,11 +212,6 @@ contains
         do ii = 1, this%nmodes
           ncols = parse_block_cols(blk, ii-1)
 
-          ! Sanity check
-          if (ncols /= this%dim + 2) then
-            call messages_input_error(namespace, 'PhotonModes', 'Incorrect number of columns')
-          end if
-
           ! Read line
           call parse_block_float(blk, ii-1, 0, this%omega(ii), units_inp%energy)  ! frequency
           call parse_block_float(blk, ii-1, 1, this%lambda(ii), units_inp%energy) ! coupling strength
@@ -214,6 +227,11 @@ contains
             call parse_block_float(blk, ii-1, 6, this%pt_momen_p0(ii))   !row, column
           end if
 
+          ! Sanity check
+          if ((ncols /= this%dim + 2) .and. (ncols /= this%dim + 4)) then
+            call messages_input_error(namespace, 'PhotonModes', 'Incorrect number of columns')
+          end if
+
           ! Normalize polarization vector
           this%pol = this%pol/sqrt(dot_product(this%pol(ii,:), this%pol(ii,:)))