controlfunction.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module controlfunction_oct_m
  use debug_oct_m
  use ext_partner_list_oct_m
  use filter_oct_m
  use global_oct_m
  use interaction_partner_oct_m
  use io_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_adv_oct_m
  use lasers_oct_m
  use math_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use tdfunction_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
 
  implicit none
 
  private
  public ::                            &
    controlfunction_t,                 &
    controlfunction_mod_init,          &
    controlfunction_mod_close,         &
    controlfunction_init,              &
    controlfunction_representation,    &
    controlfunction_mode,              &
    controlfunction_set,               &
    controlfunction_end,               &
    controlfunction_copy,              &
    controlfunction_to_h,              &
    controlfunction_to_h_val,          &
    controlfunction_write,             &
    controlfunction_diff,              &
    controlfunction_apply_envelope,    &
    controlfunction_set_fluence,       &
    controlfunction_set_alpha,         &
    controlfunction_set_rep,           &
    controlfunction_to_realtime,       &
    controlfunction_to_basis,          &
    controlfunction_prepare_initial,   &
    controlfunction_fluence,           &
    controlfunction_j2,                &
    controlfunction_get_theta,         &
    controlfunction_set_theta,         &
    controlfunction_randomize,         &
    controlfunction_update,            &
    controlfunction_number,            &
    controlfunction_bounds,            &
    controlfunction_dof,               &
    controlfunction_w0,                &
    controlfunction_alpha,             &
    controlfunction_targetfluence,     &
    controlfunction_filter,            &
    controlfunction_gradient
 
 
  integer, public, parameter ::     &
    ctr_internal               = 1, &
    ctr_fourier_series_h       = 3, &
    ctr_zero_fourier_series_h  = 4, &
    ctr_fourier_series         = 5, &
    ctr_zero_fourier_series    = 6, &
    ctr_rt                     = 7
 
 
 
  integer, parameter, public ::         &
    controlfunction_mode_none      = 0, &
    controlfunction_mode_epsilon   = 1, &
    controlfunction_mode_f         = 2
 
  type controlfunction_common_t
    private
    integer :: representation      = 0
    real(real64)   :: omegamax            = m_zero                        
    real(real64)   :: targetfluence       = m_zero                        
    logical :: fix_initial_fluence = .false.                       
    integer :: mode                = controlfunction_mode_none     
    real(real64)   :: w0                  = m_zero                        
    integer :: no_controlfunctions = 0                             !! The number of control functions to be optimized.
    real(real64),       allocatable :: alpha(:)
    type(tdf_t), allocatable :: td_penalty(:)
  end type controlfunction_common_t
 
 
  type controlfunction_t
    private
    integer :: no_controlfunctions = 0
    integer :: dim           = 0
    integer :: dof           = 0
    type(tdf_t), allocatable :: f(:)
    real(real64), allocatable :: alpha(:)
 
    integer :: current_representation = 0
 
    real(real64)   :: w0       = m_zero
    real(real64), allocatable :: u(:, :)
    real(real64), allocatable :: utransf(:, :)
    real(real64), allocatable :: utransfi(:, :)
 
    real(real64), allocatable :: theta(:)
  end type controlfunction_t
 
  logical                        :: cf_common_initialized = .false.
  type(controlfunction_common_t) :: cf_common
 
contains
 
  subroutine controlfunction_mod_init(ext_partners, namespace, dt, max_iter, mode_fixed_fluence)
    type(partner_list_t),           intent(in)    :: ext_partners
    type(namespace_t),              intent(in)    :: namespace
    real(real64),                   intent(in)    :: dt
    integer,                        intent(in)    :: max_iter
    logical,                        intent(out)   :: mode_fixed_fluence
 
    character(len=1024) :: expression
    integer :: no_lines, steps, il, idir, ncols, ipar, irow, ierr
    real(real64)   :: octpenalty
    complex(real64)   :: pol(3)
    type(block_t) :: blk
    type(lasers_t), pointer :: lasers
 
    push_sub(controlfunction_mod_init)
 
    if (.not. cf_common_initialized) then
      cf_common_initialized = .true.
    else
      message(1) = "Internal error: Cannot call controlfunction_mod_init twice."
      call messages_fatal(1, namespace=namespace)
    end if
 
    call messages_print_with_emphasis(msg="OCT: Info about control functions", namespace=namespace)
 
    !%Variable OCTControlFunctionRepresentation
    !%Type integer
    !%Section Calculation Modes::Optimal Control
    !%Default control_fourier_series_h
    !%Description
    !% If <tt>OCTControlRepresentation = control_function_parametrized</tt>, one must
    !% specify the kind of parameters that determine the control function.
    !% If <tt>OCTControlRepresentation = control_function_real_time</tt>, then this variable
    !% is ignored, and the control function is handled directly in real time.
    !%Option control_fourier_series_h 3
    !% The control function is expanded as a full Fourier series (although it must, of
    !% course, be a real function). Then, the total fluence is fixed, and a transformation
    !% to hyperspherical coordinates is done; the parameters to optimize are the hyperspherical
    !% angles.
    !%Option control_zero_fourier_series_h 4
    !% The control function is expanded as a Fourier series, but assuming (1) that the zero
    !% frequency component is zero, and (2) the control function, integrated in time, adds
    !% up to zero (this essentially means that the sum of all the cosine coefficients is zero).
    !% Then, the total fluence is fixed, and a transformation to hyperspherical coordinates is
    !% done; the parameters to optimize are the hyperspherical angles.
    !%Option control_fourier_series 5
    !% The control function is expanded as a full Fourier series (although it must, of
    !% course, be a real function). The control parameters are the coefficients of this
    !% basis-set expansion.
    !%Option control_zero_fourier_series 6
    !% The control function is expanded as a full Fourier series (although it must, of
    !% course, be a real function). The control parameters are the coefficients of this
    !% basis-set expansion. The difference with the option <tt>control_fourier_series</tt> is that
    !% (1) that the zero-frequency component is zero, and (2) the control function, integrated
    !% in time, adds up to zero (this essentially means that the sum of all the cosine
    !% coefficients is zero).
    !%Option control_rt 7
    !% (experimental)
    !%End
    call parse_variable(namespace, 'OCTControlFunctionRepresentation', ctr_rt, cf_common%representation)
    if (.not. varinfo_valid_option('OCTControlFunctionRepresentation', cf_common%representation)) then
      call messages_input_error(namespace, 'OCTControlFunctionRepresentation')
    end if
    select case (cf_common%representation)
    case (ctr_fourier_series_h)
      write(message(1), '(a)') 'Info: The OCT control functions will be represented as a Fourier series,'
      write(message(2), '(a)') '      and then a transformation to hyperspherical coordinates will be made.'
      call messages_info(2, namespace=namespace)
    case (ctr_zero_fourier_series_h)
      write(message(1), '(a)') 'Info: The OCT control functions will be represented as a Fourier series,'
      write(message(2), '(a)') '      in which (i) the zero-frequency component is assumed to be zero,'
      write(message(3), '(a)') '      and  (ii) the sum of all the cosine coefficients are zero, so that'
      write(message(4), '(a)') '      the control function starts and ends at zero.'
      write(message(5), '(a)') '      Then, a transformation to hyperspherical coordinates will be made.'
      call messages_info(5, namespace=namespace)
    case (ctr_fourier_series)
      write(message(1), '(a)') 'Info: The OCT control functions will be represented as a Fourier series.'
      call messages_info(1, namespace=namespace)
    case (ctr_zero_fourier_series)
      write(message(1), '(a)') 'Info: The OCT control functions will be represented as a Fourier series,'
      write(message(2), '(a)') '      in which the zero-frequency component is assumed to be zero,'
      write(message(3), '(a)') '      and  (ii) the sum of all the cosine coefficients are zero, so that'
      write(message(4), '(a)') '      the control function starts and ends at zero.'
      call messages_info(4, namespace=namespace)
    case (ctr_rt)
      write(message(1), '(a)') 'Info: The OCT control functions will be represented in real time.'
      call messages_info(1, namespace=namespace)
      call messages_experimental('"OCTControlFunctionRepresentation = ctr_rt"', namespace=namespace)
    case default
      write(message(1), '(a)') 'Info: The OCT control functions will be represented in real time.'
      call messages_info(1, namespace=namespace)
      call messages_experimental('"OCTControlFunctionRepresentation = ctr_rt"', namespace=namespace)
    end select
 
    !%Variable OCTControlFunctionOmegaMax
    !%Type float
    !%Section Calculation Modes::Optimal Control
    !%Default -1.0
    !%Description
    !% The Fourier series that can be used to represent the control functions must be truncated;
    !% the truncation is given by a cut-off frequency which is determined by this variable.
    !%End
    call parse_variable(namespace, 'OCTControlFunctionOmegaMax', -m_one, cf_common%omegamax)
    if (cf_common%representation /= ctr_rt) then
      write(message(1), '(a)')         'Info: The representation of the OCT control parameters will be restricted'
      write(message(2), '(a,f10.5,a)') '      with an energy cut-off of ', &
        units_from_atomic(units_out%energy, cf_common%omegamax), ' ['//trim(units_abbrev(units_out%energy)) // ']'
      call messages_info(2, namespace=namespace)
    end if
 
    !%Variable OCTFixFluenceTo
    !%Type float
    !%Section Calculation Modes::Optimal Control
    !%Default 0.0
    !%Description
    !% The algorithm tries to obtain the specified fluence for the laser field.
    !% This works only in conjunction with either the WG05 or the straight iteration scheme.
    !%
    !% If this variable is not present in the input file, by default the code will not
    !% attempt a fixed-fluence QOCT run. The same holds if the value given to this
    !% variable is exactly zero.
    !%
    !% If this variable is given a negative value, then the target fluence will be that of
    !% the initial laser pulse given as guess in the input file. Note, however, that
    !% first the code applies the envelope provided by the <tt>OCTLaserEnvelope</tt> input
    !% option, and afterwards it calculates the fluence.
    !%End
    call parse_variable(namespace, 'OCTFixFluenceTo', m_zero, cf_common%targetfluence)
 
    !%Variable OCTFixInitialFluence
    !%Type logical
    !%Section Calculation Modes::Optimal Control
    !%Default yes
    !%Description
    !% By default, when asking for a fixed-fluence optimization (<tt>OCTFixFluenceTo = whatever</tt>),
    !% the initial laser guess provided in the input file is scaled to match this
    !% fluence. However, you can force the program to use that initial laser as the initial
    !% guess, no matter the fluence, by setting <tt>OCTFixInitialFluence = no</tt>.
    !%End
    call parse_variable(namespace, 'OCTFixInitialFluence', .true., cf_common%fix_initial_fluence)
 
 
    !%Variable OCTControlFunctionType
    !%Type integer
    !%Section Calculation Modes::Optimal Control
    !%Default controlfunction_mode_epsilon
    !%Description
    !% The control function may fully determine the time-dependent form of the
    !% external field, or only the envelope function of this external field, or its phase.
    !% Or, we may have two different control functions, one of them providing the phase
    !% and the other one, the envelope.
    !%
    !% Note that, if <tt>OCTControlRepresentation = control_function_real_time</tt>, then the control
    !% function must <b>always</b> determine the full external field (THIS NEEDS TO BE FIXED).
    !%Option controlfunction_mode_epsilon   1
    !% In this case, the control function determines the full control function: namely,
    !% if we are considering the electric field of a laser, the time-dependent electric field.
    !%Option controlfunction_mode_f         2
    !% The optimization process attempts to find the best possible envelope. The full
    !% control field is this envelope times a cosine function with a "carrier" frequency.
    !% This carrier frequency is given by the carrier frequency of the <tt>TDExternalFields</tt>
    !% in the <tt>inp</tt> file.
    !%End
    call parse_variable(namespace, 'OCTControlFunctionType', controlfunction_mode_epsilon, cf_common%mode)
    if (.not. varinfo_valid_option('OCTControlFunctionType', cf_common%mode)) then
      call messages_input_error(namespace, 'OCTControlFunctionType')
    end if
    if (cf_common%representation == ctr_rt .and. (cf_common%mode /= controlfunction_mode_epsilon)) then
      call messages_input_error(namespace, 'OCTControlFunctionType')
    end if
    call messages_print_var_option('OCTControlFunctionType', cf_common%mode, namespace=namespace)
 
    lasers => list_get_lasers(ext_partners)
    if(.not. associated(lasers)) then
      message(1) = "No lasers found. Cannot proceed."
      call messages_fatal(1, namespace=namespace)
    end if
 
    ! Check that there are no complex polarization vectors.
    do il = 1, lasers%no_lasers
      pol(:) = laser_polarization(lasers%lasers(il))
      do idir = 1, 3
        if (aimag(pol(idir))**2 > 1.0e-20_real64) then
          write(message(1), '(a)') 'In QOCT runs, the polarization vector cannot be complex. Complex'
          write(message(2), '(a)') 'polarization vectors are only truly necessary if one wants a'
          write(message(3), '(a)') 'circularly / elliptically polarized laser. This concepts assumes'
          write(message(4), '(a)') 'the existence of a well defined carrier frequency (otherwise it'
          write(message(5), '(a)') 'would not make sense to speak of a fixed phase difference). So in'
          write(message(6), '(a)') 'QOCT runs it would only make sense for envelope-only optimizations.'
          write(message(7), '(a)') 'This possibility should be implemented in the future.'
          call messages_fatal(7, namespace=namespace)
        end if
      end do
    end do
 
 
    ! The laser field is defined by "td functions", as implemented in module "tdfunction_m". At this point, they
    ! can be in "non-numerical" representation (i.e. described with a set of parameters, e.g. frequency,
    ! width, etc). We need them to be in numerical form (i.e. time grid, values at the time grid).
    ! Here we do the transformation.
    ! It cannot be done before calling controlfunction_mod_init because we need to pass the omegamax value.
    do il = 1, lasers%no_lasers
      select case (cf_common%mode)
      case (controlfunction_mode_epsilon)
        call laser_to_numerical_all(lasers%lasers(il), dt, max_iter, cf_common%omegamax)
      case default
        call laser_to_numerical(lasers%lasers(il), dt, max_iter, cf_common%omegamax)
      end select
    end do
 
    ! Fix the carrier frequency
    call messages_obsolete_variable(namespace, 'OCTCarrierFrequency')
    cf_common%w0 = laser_carrier_frequency(lasers%lasers(1))
 
    ! Fix the number of control functions: if we have "traditional" QOCT (i.e. the control functions
    ! are represented directly in real time, then the number of control functions can be larger than
    ! one; it will be the number of lasers found in the input file. Otherwise, if the control function(s)
    ! are parametrized ("OCTControlRepresentation = control_function_parametrized"), we only have one
    ! control function. If there is more than one laser field in the input file, the program stops.
    if (lasers%no_lasers > 1) then
      write(message(1), '(a)') 'Currently octopus only accepts one control field.'
      call messages_fatal(1, namespace=namespace)
    end if
    cf_common%no_controlfunctions = 1
 
    mode_fixed_fluence = .false.
    select case (cf_common%representation)
    case (ctr_fourier_series_h, ctr_zero_fourier_series_h)
      if (abs(cf_common%targetfluence) <= m_epsilon) then
        write(message(1), '(a)') 'If you set "OCTControlFunctionRepresentation" to either'
        write(message(2), '(a)') '"control_fourier_series_h", or "control_zero_fourier_series_h", then the run'
        write(message(3), '(a)') 'must be done in fixed fluence mode.'
        call messages_fatal(3, namespace=namespace)
      end if
      mode_fixed_fluence = .true.
    case (ctr_fourier_series, ctr_zero_fourier_series)
      if (abs(cf_common%targetfluence) > m_epsilon) then
        write(message(1), '(a)') 'If you set "OCTControlFunctionRepresentation" to "control_fourier_series",'
        write(message(2), '(a)') 'then you cannot run in fixed fluence mode.'
        call messages_fatal(2, namespace=namespace)
      end if
      mode_fixed_fluence = .false.
    case default
      if (abs(cf_common%targetfluence) > m_epsilon) mode_fixed_fluence = .true.
    end select
 
 
    !%Variable OCTPenalty
    !%Type float
    !%Section Calculation Modes::Optimal Control
    !%Default 1.0
    !%Description
    !% The variable specifies the value of the penalty factor for the
    !% integrated field strength (fluence). Large value = small fluence.
    !% A transient shape can be specified using the block <tt>OCTLaserEnvelope</tt>.
    !% In this case <tt>OCTPenalty</tt> is multiplied with time-dependent function.
    !% The value depends on the coupling between the states. A good start might be a
    !% value from 0.1 (strong fields) to 10 (weak fields).
    !%
    !% Note that if there are several control functions, one can specify this
    !% variable as a one-line code, each column being the penalty factor for each
    !% of the control functions. Make sure that the number of columns is equal to the
    !% number of control functions. If it is not a block, all control functions will
    !% have the same penalty factor.
    !%
    !% All penalty factors must be positive.
    !%End
    safe_allocate(cf_common%alpha(1:cf_common%no_controlfunctions))
    cf_common%alpha = m_zero
    if (parse_block(namespace, 'OCTPenalty', blk) == 0) then
      ! We have a block
      ncols = parse_block_cols(blk, 0)
      if (ncols /= cf_common%no_controlfunctions) then
        call messages_input_error(namespace, 'OCTPenalty')
      else
        do ipar = 1, ncols
          call parse_block_float(blk, 0, ipar - 1, cf_common%alpha(ipar))
          if (cf_common%alpha(ipar) <= m_zero) call messages_input_error(namespace, 'OCTPenalty')
        end do
      end if
    else
      ! We have the same penalty for all the control functions.
      call parse_variable(namespace, 'OCTPenalty', m_one, octpenalty)
      cf_common%alpha(1:cf_common%no_controlfunctions) = octpenalty
    end if
 
 
    !%Variable OCTLaserEnvelope
    !%Type block
    !%Section Calculation Modes::Optimal Control
    !%Description
    !% Often a pre-defined time-dependent envelope on the control function is desired.
    !% This can be achieved by making the penalty factor time-dependent.
    !% Here, you may specify the required time-dependent envelope.
    !%
    !% It is possible to choose different envelopes for different control functions.
    !% There should be one line for each control function. Each line should
    !% have only one element: a string with the name of a time-dependent function,
    !% that should be correspondingly defined in a <tt>TDFunctions</tt> block.
    !%End
    steps = max_iter
    safe_allocate(cf_common%td_penalty(1:cf_common%no_controlfunctions))
 
    if (parse_block(namespace, 'OCTLaserEnvelope', blk) == 0) then
 
      ! Cannot have this unless we have the "usual" controlfunction_mode_epsilon.
      if (cf_common%mode /= controlfunction_mode_epsilon) then
        write(message(1),'(a)') 'The block "OCTLaserEnvelope" is only compatible with the option'
        write(message(2),'(a)') '"OCTControlFunctionType = controlfunction_mode_epsilon".'
        call messages_fatal(2, namespace=namespace)
      end if
 
      no_lines = parse_block_n(blk)
      if (no_lines /= cf_common%no_controlfunctions) call messages_input_error(namespace, 'OCTLaserEnvelope')
 
      do irow = 1, no_lines
        call parse_block_string(blk, irow - 1, 0, expression)
        call parse_block_end(blk)
        call tdf_read(cf_common%td_penalty(irow), namespace, trim(expression), ierr)
        if (ierr.ne.0) then
          message(1) = 'Time-dependent function "'//trim(expression)//'" could not be read from inp file.'
          call messages_fatal(1, namespace=namespace)
        end if
        call tdf_to_numerical(cf_common%td_penalty(irow), steps, dt, cf_common%omegamax)
        ierr = parse_block(namespace, 'OCTLaserEnvelope', blk)
      end do
 
      call parse_block_end(blk)
    else
      do ipar = 1, cf_common%no_controlfunctions
        call tdf_init_numerical(cf_common%td_penalty(ipar), steps, dt, -m_one, initval = m_one)
      end do
    end if
 
    call messages_print_with_emphasis(namespace=namespace)
 
    pop_sub(controlfunction_mod_init)
  end subroutine controlfunction_mod_init
 
  ! ---------------------------------------------------------
  subroutine controlfunction_init(cp, dt, ntiter)
    type(controlfunction_t), intent(inout) :: cp
    real(real64),            intent(in)    :: dt
    integer,                 intent(in)    :: ntiter
 
    integer :: ipar
 
    push_sub(controlfunction_init)
 
    cp%w0                  = cf_common%w0
    cp%no_controlfunctions = cf_common%no_controlfunctions
    cp%current_representation = ctr_internal
    safe_allocate_source_a(cp%alpha, cf_common%alpha)
 
    safe_allocate(cp%f(1:cp%no_controlfunctions))
    do ipar = 1, cp%no_controlfunctions
      call tdf_init_numerical(cp%f(ipar), ntiter, dt, cf_common%omegamax)
    end do
 
    ! If the control function is represented directly in real time, the "dimension" (cp%dim) is
    ! the number of values that represent the function on the discretized time-axis.
    !
    ! If the control function is parametrized, up to now (in the future this might change), all
    ! parametrizations are based on a previous basis-set expansion (sine-Fourier series, or "normal"
    ! Fourier series with or without the zero term). For the representations whose name ends in "_h",
    ! the parameters are not directly the coefficients of the control function in this basis-set
    ! expansion, but are constructed from them (e.g. by performing a coordinate transformation to
    ! hyperspherical coordinates). The "dimension" (cp%dim) is the dimension of this basis set.
    select case (cf_common%representation)
    case (ctr_internal)
      cp%dim = ntiter + 1
    case (ctr_fourier_series_h, ctr_fourier_series)
      ! If nf is the number of frequencies, we will have nf-1 non-zero "sines", nf-1 non-zero "cosines",
      ! and the zero-frequency component. Total, 2*(nf-1)+1
      cp%dim = 2 * (tdf_nfreqs(cp%f(1)) - 1) + 1
    case (ctr_zero_fourier_series, ctr_zero_fourier_series_h)
      ! If nf is the number of frequencies, we will have nf-1 non-zero "sines", nf-1 non-zero "cosines",
      ! but no zero-frequency component. Total, 2*(nf-1)
      cp%dim = 2 * (tdf_nfreqs(cp%f(1)) - 1)
    case (ctr_rt)
      cp%dim = ntiter + 1
    case default
      message(1) = "Internal error: invalid representation."
      call messages_fatal(1)
    end select
 
    ! The "degrees of freedom" cp%dof is the number of parameters that define the control function.
    ! (if it is represented directly in real time, this would be meaningless, but we put the number of
    ! control functions, times the "dimension", which in this case is the number of time discretization
    ! points). This is not equal to the dimension of the basis set employed (cp%dim), because we may
    ! add further constraints, and do a coordinate transformation to account for them.
    select case (cf_common%representation)
    case (ctr_internal)
      cp%dof = cp%no_controlfunctions * cp%dim
    case (ctr_fourier_series_h)
      ! The number of degrees of freedom is one fewer than the number of basis coefficients, since we
      ! add the constraint of fixed fluence.
      cp%dof = cp%dim - 1
    case (ctr_zero_fourier_series_h)
      ! The number of degrees of freedom is one fewer than the number of basis coefficients, since we
      ! add (1) the constraint of fixed fluence, and (2) the constraint of the field starting and
      ! ending at zero, which amounts to having all the cosine coefficients summing up to zero.
      cp%dof = cp%dim - 2
    case (ctr_fourier_series)
      ! In this case, we have no constraints: the dof is equal to the dimension of the basis set, since
      ! the parameters are directly the coefficients of the basis-set expansion.
      cp%dof = cp%dim
    case (ctr_zero_fourier_series)
      ! The number of degrees of freedom is reduced by one, since we add the constraint forcing the
      ! the field to start and end at zero, which amounts to having all the cosine coefficients
      ! summing up to zero.
      cp%dof = cp%dim - 1
    case (ctr_rt)
      cp%dof = cp%no_controlfunctions * cp%dim
    end select
 
    if (cp%dof <= 0) then
      write(message(1),'(a)') 'The number of degrees of freedom used to describe the control function'
      write(message(2),'(a)') 'is less than or equal to zero. This should not happen. Please review your input file.'
      call messages_fatal(2)
    else
      if (cf_common%representation /= ctr_internal) then
        write(message(1), '(a,i6,a)') 'The parametrization of the control functions makes use of ', cp%dim, ' basis'
        write(message(2), '(a,i6,a)') 'functions and ', cp%dof, ' degrees of freedom.'
        call messages_info(2)
        safe_allocate(cp%theta(1:cp%dof))
        cp%theta = m_zero
      end if
    end if
 
    ! Construct the transformation matrix, if needed.
    call controlfunction_trans_matrix(cp)
 
    pop_sub(controlfunction_init)
  end subroutine controlfunction_init
  ! ---------------------------------------------------------
 
 
 
  subroutine controlfunction_set(cp, ext_partners)
    type(controlfunction_t), intent(inout) :: cp
    type(partner_list_t),       intent(in) :: ext_partners
 
    integer :: ipar
    type(lasers_t), pointer :: lasers
 
    push_sub(controlfunction_set)
 
    lasers => list_get_lasers(ext_partners)
    assert(associated(lasers))
 
    select case (cf_common%mode)
    case (controlfunction_mode_epsilon, controlfunction_mode_f)
      do ipar = 1, cp%no_controlfunctions
        call tdf_end(cp%f(ipar))
        call laser_get_f(lasers%lasers(ipar), cp%f(ipar))
      end do
    end select
 
    pop_sub(controlfunction_set)
  end subroutine controlfunction_set
  ! ---------------------------------------------------------
 
 
  integer pure function controlfunction_representation()
    controlfunction_representation = cf_common%representation
  end function controlfunction_representation
  ! ---------------------------------------------------------
 
 
  integer pure function controlfunction_mode()
    controlfunction_mode = cf_common%mode
  end function controlfunction_mode
  ! ---------------------------------------------------------
 
 
  subroutine controlfunction_prepare_initial(par)
    type(controlfunction_t), intent(inout) :: par
 
    push_sub(controlfunction_prepare_initial)
 
    call controlfunction_apply_envelope(par)
 
    if (abs(cf_common%targetfluence) > m_epsilon) then
      if (cf_common%targetfluence < m_zero) then
        cf_common%targetfluence = controlfunction_fluence(par)
        write(message(1), '(a)')         'Info: The QOCT run will attempt to find a solution with the same'
        write(message(2), '(a,f10.5,a)') '      fluence as the input external fields: F = ', &
          cf_common%targetfluence, ' a.u.'
      else
        write(message(1), '(a)')         'Info: The QOCT run will attempt to find a solution with a predefined'
        write(message(2), '(a,f10.5,a)') '      fluence: F = ', cf_common%targetfluence, ' a.u.'
      end if
      call messages_info(2)
      if (cf_common%fix_initial_fluence) call controlfunction_set_fluence(par)
    end if
 
    ! Move to the "native" representation, if necessary.
    call controlfunction_set_rep(par)
 
    pop_sub(controlfunction_prepare_initial)
  end subroutine controlfunction_prepare_initial
  ! ---------------------------------------------------------
 
 
  subroutine controlfunction_set_rep(par)
    type(controlfunction_t), intent(inout) :: par
 
    push_sub(controlfunction_set_rep)
 
    if (par%current_representation /= cf_common%representation) then
      if (cf_common%representation == ctr_internal) then
        call controlfunction_to_realtime(par)
      else
        call controlfunction_to_basis(par)
      end if
    end if
 
    pop_sub(controlfunction_set_rep)
  end subroutine controlfunction_set_rep
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine controlfunction_to_basis(par)
    type(controlfunction_t), intent(inout) :: par
 
    integer :: ipar
 
    push_sub(controlfunction_to_basis)
 
    if (par%current_representation == ctr_internal) then
      select case (cf_common%representation)
      case (ctr_rt)
        par%current_representation = ctr_rt
        call controlfunction_basis_to_theta(par)
      case (ctr_fourier_series_h)
        do ipar = 1, par%no_controlfunctions
          call tdf_numerical_to_fourier(par%f(ipar))
        end do
        par%current_representation = ctr_fourier_series_h
        call controlfunction_basis_to_theta(par)
      case (ctr_zero_fourier_series_h)
        do ipar = 1, par%no_controlfunctions
          call tdf_numerical_to_zerofourier(par%f(ipar))
        end do
        par%current_representation = ctr_zero_fourier_series_h
        call controlfunction_basis_to_theta(par)
      case (ctr_fourier_series)
        do ipar = 1, par%no_controlfunctions
          call tdf_numerical_to_fourier(par%f(ipar))
        end do
        par%current_representation = ctr_fourier_series
        call controlfunction_basis_to_theta(par)
      case (ctr_zero_fourier_series)
        do ipar = 1, par%no_controlfunctions
          call tdf_numerical_to_zerofourier(par%f(ipar))
        end do
        par%current_representation = ctr_zero_fourier_series
        call controlfunction_basis_to_theta(par)
      end select
    end if
 
    pop_sub(controlfunction_to_basis)
  end subroutine controlfunction_to_basis
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine controlfunction_to_realtime(par)
    type(controlfunction_t), intent(inout) :: par
 
    integer :: ipar
 
    push_sub(controlfunction_to_realtime)
 
    select case (par%current_representation)
    case (ctr_internal)
      pop_sub(controlfunction_to_realtime)
      return
    case (ctr_fourier_series_h)
      call controlfunction_theta_to_basis(par)
      do ipar = 1, par%no_controlfunctions
        call tdf_fourier_to_numerical(par%f(ipar))
      end do
    case (ctr_zero_fourier_series_h)
      call controlfunction_theta_to_basis(par)
      do ipar = 1, par%no_controlfunctions
        call tdf_zerofourier_to_numerical(par%f(ipar))
      end do
    case (ctr_fourier_series)
      call controlfunction_theta_to_basis(par)
      do ipar = 1, par%no_controlfunctions
        call tdf_fourier_to_numerical(par%f(ipar))
      end do
    case (ctr_zero_fourier_series)
      call controlfunction_theta_to_basis(par)
      do ipar = 1, par%no_controlfunctions
        call tdf_zerofourier_to_numerical(par%f(ipar))
      end do
    case (ctr_rt)
      call controlfunction_theta_to_basis(par)
    end select
 
    par%current_representation = ctr_internal
    pop_sub(controlfunction_to_realtime)
  end subroutine controlfunction_to_realtime
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  real(real64) function controlfunction_diff(pp, qq) result(res)
    type(controlfunction_t), intent(in) :: pp, qq
 
    integer :: ipar
 
    push_sub(controlfunction_diff)
 
    assert(pp%current_representation == qq%current_representation)
 
    res = m_zero
    do ipar = 1, pp%no_controlfunctions
      res = res + tdf_diff(pp%f(ipar), qq%f(ipar))
    end do
 
    pop_sub(controlfunction_diff)
  end function controlfunction_diff
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine controlfunction_apply_envelope(cp)
    type(controlfunction_t), intent(inout) :: cp
    integer :: ipar, iter
 
    push_sub(controlfunction_apply_envelope)
 
    ! Do not apply the envelope if the control functions are represented as a sine-Fourier series.
    if (cf_common%representation == ctr_rt) then
      call controlfunction_to_realtime(cp)
      do ipar = 1, cp%no_controlfunctions
        do iter = 1, tdf_niter(cp%f(ipar)) + 1
          call tdf_set_numerical(cp%f(ipar), iter, tdf(cp%f(ipar), iter) / tdf(cf_common%td_penalty(ipar), iter))
        end do
      end do
      call controlfunction_to_basis(cp)
    end if
 
    pop_sub(controlfunction_apply_envelope)
  end subroutine controlfunction_apply_envelope
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine controlfunction_to_h(cp, ext_partners)
    type(controlfunction_t), intent(in) :: cp
    type(partner_list_t),    intent(in) :: ext_partners
 
    integer :: ipar
    type(controlfunction_t) :: par
    type(lasers_t), pointer :: lasers
 
    push_sub(controlfunction_to_h)
 
    call controlfunction_copy(par, cp)
    call controlfunction_to_realtime(par)
 
    lasers => list_get_lasers(ext_partners)
    assert(associated(lasers))
 
    select case (cf_common%mode)
    case (controlfunction_mode_epsilon, controlfunction_mode_f)
      do ipar = 1, cp%no_controlfunctions
        call laser_set_f(lasers%lasers(ipar), par%f(ipar))
      end do
    end select
 
    call controlfunction_end(par)
 
    pop_sub(controlfunction_to_h)
  end subroutine controlfunction_to_h
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine controlfunction_to_h_val(cp, ext_partners, val)
    type(controlfunction_t), intent(in) :: cp
    type(partner_list_t),    intent(in) :: ext_partners
    integer,                 intent(in) :: val
 
    integer :: ipar
    type(lasers_t), pointer :: lasers
 
    push_sub(controlfunction_to_h_val)
 
    lasers => list_get_lasers(ext_partners)
    assert(associated(lasers))
 
    do ipar = 1, cp%no_controlfunctions
      call laser_set_f_value(lasers%lasers(ipar), val, tdf(cp%f(ipar), val))
    end do
 
    pop_sub(controlfunction_to_h_val)
  end subroutine controlfunction_to_h_val
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine controlfunction_end(cp)
    type(controlfunction_t), intent(inout) :: cp
    integer :: ipar
 
    push_sub(controlfunction_end)
 
    if (allocated(cp%f)) then
      do ipar = 1, cp%no_controlfunctions
        call tdf_end(cp%f(ipar))
      end do
    end if
    safe_deallocate_a(cp%f)
    safe_deallocate_a(cp%alpha)
    safe_deallocate_a(cp%u)
    safe_deallocate_a(cp%utransf)
    safe_deallocate_a(cp%utransfi)
    safe_deallocate_a(cp%theta)
 
    pop_sub(controlfunction_end)
  end subroutine controlfunction_end
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine controlfunction_write(filename, cp, namespace)
    character(len=*), intent(in) :: filename
    type(controlfunction_t), intent(in) :: cp
    type(namespace_t),       intent(in) :: namespace
 
    integer :: iter, ipar, ifreq, iunit, niter, nfreqs, idof
    real(real64) :: time, wmax, dw, ww, wa, wb, dt
    real(real64), allocatable :: func(:, :)
    complex(real64) :: ft, ez, ezdt
    character(len=2) :: digit
    type(controlfunction_t) :: par
 
    if (.not. mpi_grp_is_root(mpi_world)) return
 
    push_sub(controlfunction_write)
 
    call io_mkdir(trim(filename), namespace)
 
    call controlfunction_copy(par, cp)
    call controlfunction_to_realtime(par)
 
    iunit = io_open(trim(filename)//'/Fluence', namespace, action='write')
    write(iunit, '(a,es20.8e3)') 'Fluence = ', controlfunction_fluence(par)
    call io_close(iunit)
 
    niter = tdf_niter(par%f(1))
    safe_allocate(func(1:niter + 1, 1:cp%no_controlfunctions))
 
    select case (cf_common%mode)
    case (controlfunction_mode_epsilon)
 
      do ipar = 1, cp%no_controlfunctions
        if (cp%no_controlfunctions > 1) then
          write(digit,'(i2.2)') ipar
          iunit = io_open(trim(filename)//'/cp-'//digit, namespace, action='write')
        else
          iunit = io_open(trim(filename)//'/cp', namespace, action='write')
        end if
        write(iunit,'(2a20)') '#       t [a.u]      ', '        e(t)         '
        do iter = 1, tdf_niter(par%f(ipar)) + 1
          time = (iter - 1) * tdf_dt(par%f(ipar))
          write(iunit, '(2es20.8e3)') time, tdf(par%f(ipar), iter)
          func(iter, ipar) = tdf(par%f(ipar), time)
        end do
        call io_close(iunit)
      end do
 
    case (controlfunction_mode_f)
 
      do ipar = 1, cp%no_controlfunctions
        if (cp%no_controlfunctions > 1) then
          write(digit,'(i2.2)') ipar
          iunit = io_open(trim(filename)//'/cp-'//digit, namespace, action='write')
        else
          iunit = io_open(trim(filename)//'/cp', namespace, action='write')
        end if
        write(iunit,'(3a20)') '#       t [a.u]      ', '        e(t)         ', '        f(t)         '
        do iter = 1, tdf_niter(par%f(ipar)) + 1
          time = (iter - 1) * tdf_dt(par%f(ipar))
          write(iunit, '(3es20.8e3)') time, tdf(par%f(ipar), time) * cos(par%w0 * time), tdf(par%f(ipar), time)
          func(iter, ipar) = tdf(par%f(ipar), time) * cos(par%w0 * time)
        end do
        call io_close(iunit)
      end do
 
    end select
 
 
    !Now, the Fourier transforms.
    select case (cf_common%mode)
    case (controlfunction_mode_epsilon)
 
      do ipar = 1, cp%no_controlfunctions
        if (cp%no_controlfunctions > 1) then
          write(digit,'(i2.2)') ipar
          iunit = io_open(trim(filename)//'/cpw-'//digit, namespace, action='write')
        else
          iunit = io_open(trim(filename)//'/cpw', namespace, action='write')
        end if
        write(iunit,'(3a20)') '#       w [a.u]      ', '      Re[e(w)]       ', &
          '      Im[e(w)]       '
 
        nfreqs = 1000
        wa = m_zero
        wb = m_three ! hard-coded to three atomic units... this should be improved.
        wmax = wb
        dw = wmax / (nfreqs - 1)
        dt = tdf_dt(par%f(1))
 
        do ifreq = 1, nfreqs
          ww = wa + (ifreq - 1) * dw
          ft = m_z0
          ez = m_z1
          ezdt = exp(m_zi * ww * tdf_dt(par%f(ipar)))
          do iter = 1, niter + 1
            time = (iter - 1) * dt
            ft = ft + func(iter, ipar) * ez
            ez = ez * ezdt
          end do
          ft = ft * dt
          write(iunit,'(3es20.8e3)') ww, real(ft, real64) , aimag(ft)
        end do
        call io_close(iunit)
      end do
 
 
    case (controlfunction_mode_f)
      iunit = io_open(trim(filename)//'/cpw', namespace, action='write')
      write(iunit,'(3a20)') '#       w [a.u]      ', '      Re[e(w)]       ', &
        '      Im[e(w)]       '
 
      nfreqs = 1000
      wa = cp%w0 - m_three * cf_common%omegamax
      wb = cp%w0 + m_three * cf_common%omegamax
      wmax = 6.0_real64*cf_common%omegamax
      dw = wmax/(nfreqs-1)
      dt = tdf_dt(par%f(1))
 
      do ifreq = 1, nfreqs
        ww = wa + (ifreq - 1) * dw
        ft = m_z0
        ez = m_z1
        ezdt = exp(m_zi * ww * tdf_dt(par%f(1)))
        do iter = 1, niter + 1
          time = (iter - 1) * dt
          ft = ft + func(iter, 1) * ez
          ez = ez * ezdt
        end do
        ft = ft * dt
        write(iunit,'(3es20.8e3)') ww, real(ft, real64) , aimag(ft)
      end do
 
      call io_close(iunit)
    end select
 
    ! Now, in case of a parametrized control function, the parameters.
    if (cf_common%representation /= ctr_rt) then
      iunit = io_open(trim(filename)//'/theta', namespace, action='write')
      do idof = 1, par%dof
        write(iunit,'(i5,es20.8e3)') idof, par%theta(idof)
      end do
      call io_close(iunit)
    end if
 
    call controlfunction_end(par)
    pop_sub(controlfunction_write)
  end subroutine controlfunction_write
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  ! Gets the fluence of the laser field, defined as:
  ! controlfunction_fluence = \sum_i^{no_controlfunctions} \integrate_0^T |epsilon(t)|^2
  ! ---------------------------------------------------------
  real(real64) function controlfunction_fluence(par)
    type(controlfunction_t), intent(in) :: par
    type(controlfunction_t)             :: par_
    integer :: ipar
    type(tdf_t) :: ff
    push_sub(controlfunction_fluence)
 
    call controlfunction_copy(par_, par)
    call controlfunction_to_realtime(par_)
 
    controlfunction_fluence = m_zero
 
    select case (cf_common%mode)
    case (controlfunction_mode_epsilon)
      do ipar = 1, par_%no_controlfunctions
        controlfunction_fluence = controlfunction_fluence + tdf_dot_product(par_%f(ipar), par_%f(ipar))
      end do
    case (controlfunction_mode_f)
      do ipar = 1, par%no_controlfunctions
        call tdf_init(ff)
        call tdf_copy(ff, par_%f(ipar))
        call tdf_cosine_multiply(par%w0, ff)
        controlfunction_fluence = controlfunction_fluence + tdf_dot_product(ff, ff)
        call tdf_end(ff)
      end do
    end select
 
    call controlfunction_end(par_)
    pop_sub(controlfunction_fluence)
  end function controlfunction_fluence
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  ! Gets the J2 functional (which is the fluence, but weighted
  ! by a penalty function.
  ! ---------------------------------------------------------
  real(real64) function controlfunction_j2(par) result(j2)
    type(controlfunction_t), intent(in) :: par
    type(controlfunction_t)             :: par_
    integer :: iter, ipar
    real(real64)   :: time, integral, fi, tdp
    type(tdf_t) :: ff
 
    push_sub(controlfunction_j2)
 
    assert(par%current_representation == cf_common%representation)
 
    call controlfunction_copy(par_, par)
    call controlfunction_to_realtime(par_)
 
    integral = m_zero
    select case (cf_common%mode)
    case (controlfunction_mode_epsilon)
      do ipar = 1, par_%no_controlfunctions
        call tdf_init(ff)
        call tdf_copy(ff, par_%f(ipar))
        do iter = 1, tdf_niter(ff) + 1
          time = (iter - 1) * tdf_dt(ff)
          fi = tdf(par_%f(ipar), iter)
          tdp = sqrt(real(tdf(cf_common%td_penalty(ipar), iter), real64))
          call tdf_set_numerical(ff, iter, fi * tdp)
        end do
        integral = integral + tdf_dot_product(ff, ff)
        call tdf_end(ff)
      end do
    case (controlfunction_mode_f)
      do ipar = 1, par_%no_controlfunctions
        call tdf_init(ff)
        call tdf_copy(ff, par_%f(ipar))
        do iter = 1, tdf_niter(ff) + 1
          time = (iter - 1) * tdf_dt(ff)
          fi = tdf(par_%f(ipar), iter)
          tdp = sqrt(real(tdf(cf_common%td_penalty(ipar), iter), real64))
          call tdf_set_numerical(ff, iter, fi * tdp * cos(par_%w0 * time))
        end do
        integral = integral + tdf_dot_product(ff, ff)
        call tdf_end(ff)
      end do
    end select
 
    j2 = - par_%alpha(1) * (integral - cf_common%targetfluence)
 
    call controlfunction_end(par_)
    pop_sub(controlfunction_j2)
  end function controlfunction_j2
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine controlfunction_set_fluence(par)
    type(controlfunction_t), intent(inout) :: par
    real(real64)   :: old_fluence
    integer :: ipar
 
    push_sub(controlfunction_set_fluence)
 
    call controlfunction_to_realtime(par)
    old_fluence = controlfunction_fluence(par)
    do ipar = 1, par%no_controlfunctions
      call tdf_scalar_multiply( sqrt(cf_common%targetfluence / old_fluence), par%f(ipar))
    end do
    call controlfunction_to_basis(par)
 
    pop_sub(controlfunction_set_fluence)
  end subroutine controlfunction_set_fluence
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine controlfunction_set_alpha(par, alpha)
    type(controlfunction_t), intent(inout) :: par
 
    real(real64), intent(in) :: alpha
 
    push_sub(controlfunction_set_alpha)
 
    par%alpha(:) = alpha
 
    pop_sub(controlfunction_set_alpha)
  end subroutine controlfunction_set_alpha
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine controlfunction_copy(cp_out, cp_in)
    type(controlfunction_t), intent(inout) :: cp_out
    type(controlfunction_t), intent(in)    :: cp_in
 
    integer :: ipar
 
    push_sub(controlfunction_copy)
 
    call controlfunction_end(cp_out)
 
    cp_out%no_controlfunctions = cp_in%no_controlfunctions
    cp_out%dim = cp_in%dim
    cp_out%dof = cp_in%dof
    cp_out%current_representation = cp_in%current_representation
    cp_out%w0 = cp_in%w0
 
    safe_allocate_source_a(cp_out%alpha, cp_in%alpha)
    safe_allocate(cp_out%f(1:cp_out%no_controlfunctions))
 
    do ipar = 1, cp_in%no_controlfunctions
      call tdf_init(cp_out%f(ipar))
      call tdf_copy(cp_out%f(ipar), cp_in%f(ipar))
    end do
 
    safe_allocate_source_a(cp_out%u, cp_in%u)
    safe_allocate_source_a(cp_out%utransf, cp_in%utransf)
    safe_allocate_source_a(cp_out%utransfi, cp_in%utransfi)
    safe_allocate_source_a(cp_out%theta, cp_in%theta)
 
    pop_sub(controlfunction_copy)
  end subroutine controlfunction_copy
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine controlfunction_randomize(par)
    type(controlfunction_t), intent(inout) :: par
 
    integer :: ipar
 
    push_sub(controlfunction_randomize)
 
    assert(cf_common%representation /= ctr_internal)
 
    call controlfunction_set_rep(par)
 
    do ipar = 1, par%no_controlfunctions
      call tdf_set_random(par%f(ipar))
    end do
    call controlfunction_basis_to_theta(par)
 
    pop_sub(controlfunction_randomize)
  end subroutine controlfunction_randomize
  ! ---------------------------------------------------------
 
 
 
  subroutine controlfunction_update(cp, cpp, dir, iter, mu, dd, dq)
    type(controlfunction_t), intent(inout) :: cp
    type(controlfunction_t), intent(in)    :: cpp
    character(len=1),        intent(in)    :: dir
    integer,                 intent(in)    :: iter
    real(real64),            intent(in)    :: mu
    real(real64),            intent(in)    :: dd(:)
    complex(real64),         intent(in)    :: dq(:)
 
    real(real64) :: val
    integer :: ipar
 
    push_sub(controlfunction_update)
 
    select case (dir)
    case ('f')
      do ipar = 1, cp%no_controlfunctions
        val = dd(ipar) / ( tdf(cf_common%td_penalty(ipar), iter) - m_two * aimag(dq(ipar)))
        val = (m_one - mu) * tdf(cpp%f(ipar), iter) + mu * val
        call tdf_set_numerical(cp%f(ipar), iter, val)
        if (iter + 1 <= tdf_niter(cp%f(ipar)) + 1)  call tdf_set_numerical(cp%f(ipar), iter+1, val)
        if (iter + 2 <= tdf_niter(cp%f(ipar)) + 1)  call tdf_set_numerical(cp%f(ipar), iter+2, val)
      end do
 
    case ('b')
      do ipar = 1, cp%no_controlfunctions
        val = dd(ipar) / ( tdf(cf_common%td_penalty(ipar), iter + 1) - m_two * aimag(dq(ipar)))
        val = (m_one - mu) * tdf(cpp%f(ipar), iter + 1) + mu * val
        call tdf_set_numerical(cp%f(ipar), iter + 1, val)
        if (iter > 0) call tdf_set_numerical(cp%f(ipar), iter, val)
        if (iter - 1 > 0) call tdf_set_numerical(cp%f(ipar), iter-1, val)
      end do
    end select
 
    pop_sub(controlfunction_update)
  end subroutine controlfunction_update
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  real(real64) pure function controlfunction_alpha(par, ipar)
    type(controlfunction_t), intent(in) :: par
    integer,                 intent(in) :: ipar
    controlfunction_alpha = par%alpha(ipar)
  end function controlfunction_alpha
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  real(real64) pure function controlfunction_targetfluence()
    controlfunction_targetfluence = cf_common%targetfluence
  end function controlfunction_targetfluence
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  integer pure function controlfunction_number(par)
    type(controlfunction_t), intent(in) :: par
    controlfunction_number = par%no_controlfunctions
  end function controlfunction_number
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine controlfunction_bounds(par, lower_bounds, upper_bounds)
    type(controlfunction_t), intent(in)  :: par
    real(real64),            intent(out) :: lower_bounds(:)
    real(real64),            intent(out) :: upper_bounds(:)
 
    push_sub(controlfunction_bounds)
 
    select case (cf_common%representation)
    case (ctr_zero_fourier_series_h, ctr_fourier_series_h)
      lower_bounds(1:par%dof-1) = m_zero
      lower_bounds(par%dof) = -m_pi
      upper_bounds(1:par%dof) = m_pi
    case default
      lower_bounds = -huge(m_zero)
      lower_bounds = -huge(m_zero)
      upper_bounds =  huge(m_zero)
      upper_bounds =  huge(m_zero)
    end select
 
    pop_sub(controlfunction_bounds)
  end subroutine controlfunction_bounds
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  integer pure function controlfunction_dof(par)
    type(controlfunction_t), intent(in) :: par
    controlfunction_dof = par%dof
  end function controlfunction_dof
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  real(real64) pure function controlfunction_w0(par)
    type(controlfunction_t), intent(in) :: par
    controlfunction_w0 = par%w0
  end function controlfunction_w0
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine controlfunction_filter(par, filter)
    type(controlfunction_t), intent(inout) :: par
    type(filter_t),          intent(inout) :: filter
 
    integer :: ipar
 
    push_sub(controlfunction_filter)
 
    call controlfunction_to_realtime(par)
 
    do ipar = 1, par%no_controlfunctions
      call filter_apply(par%f(ipar), filter)
    end do
 
    call controlfunction_to_basis(par)
 
    pop_sub(controlfunction_filter)
  end subroutine controlfunction_filter
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine controlfunction_mod_close()
    integer :: ipar
 
    push_sub(controlfunction_mod_close)
 
    if (.not. cf_common_initialized) then
      message(1) = "Internal error: Cannot call controlfunction_mod_close when not initialized."
      call messages_fatal(1)
    end if
 
    cf_common_initialized = .false.
    safe_deallocate_a(cf_common%alpha)
 
    do ipar = 1, cf_common%no_controlfunctions
      call tdf_end(cf_common%td_penalty(ipar))
    end do
 
    safe_deallocate_a(cf_common%td_penalty)
 
    pop_sub(controlfunction_mod_close)
  end subroutine controlfunction_mod_close
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine controlfunction_deltaedeltau(par, dedu)
    type(controlfunction_t), intent(in)    :: par
    real(real64), intent(inout)                   :: dedu(:, :)
 
    integer :: i
    real(real64) :: rr
    real(real64), allocatable :: grad_matrix(:, :), dedv(:, :)
 
    push_sub(controlfunction_deltaedeltau)
 
    select case (cf_common%representation)
 
    case (ctr_rt)
      dedu = diagonal_matrix(par%dim, m_one)
 
    case (ctr_fourier_series_h)
      safe_allocate(grad_matrix(1:par%dim - 1, 1:par%dim))
      rr = sqrt(cf_common%targetfluence)
      call hypersphere_grad_matrix(grad_matrix, rr, par%theta)
      dedu = matmul(grad_matrix, transpose(par%utransfi))
      safe_deallocate_a(grad_matrix)
 
    case (ctr_zero_fourier_series_h)
      safe_allocate(dedv(1:par%dim-1, 1:par%dim))
      dedv = m_zero
      do i = 1, (par%dim-1)/2
        dedv(i, 1) = - m_one
        dedv(i, i+1) = m_one
      end do
      do i = (par%dim-1)/2+1, par%dim-1
        dedv(i, i+1) = m_one
      end do
      safe_allocate(grad_matrix(1:par%dim - 2, 1:par%dim - 1))
      rr = sqrt(cf_common%targetfluence)
      call hypersphere_grad_matrix(grad_matrix, rr, par%theta)
      dedu = matmul(grad_matrix, matmul(transpose(par%utransfi), dedv))
      safe_deallocate_a(grad_matrix)
      safe_deallocate_a(dedv)
 
    case (ctr_fourier_series)
      dedu = diagonal_matrix(par%dim, m_one)
 
    case (ctr_zero_fourier_series)
      dedu = m_zero
      do i = 1, par%dof/2
        dedu(i, 1) = - m_one
        dedu(i, i+1) = m_one
      end do
      do i = par%dof/2+1, par%dof
        dedu(i, i+1) = m_one
      end do
 
    end select
 
    pop_sub(controlfunction_deltaedeltau)
  end subroutine controlfunction_deltaedeltau
  ! ---------------------------------------------------------
 
 
  subroutine controlfunction_der(par, depsilon, i)
    type(controlfunction_t), intent(in)    :: par
    type(tdf_t),             intent(inout) :: depsilon
    integer,                 intent(in)    :: i
 
    real(real64), allocatable :: dedu(:, :)
    push_sub(controlfunction_der)
    assert( i > 0 .and. i <= par%dof)
 
    safe_allocate(dedu(1:par%dof, 1:par%dim))
    call controlfunction_deltaedeltau(par, dedu)
    call tdf_set_numerical(depsilon, dedu(i, 1:par%dim))
    call tdf_to_numerical(depsilon)
    if (controlfunction_mode() == controlfunction_mode_f) call tdf_cosine_multiply(par%w0, depsilon)
 
    safe_deallocate_a(dedu)
    pop_sub(controlfunction_der)
  end subroutine controlfunction_der
 
 
  subroutine controlfunction_gradient(par, par_output, grad)
    type(controlfunction_t), intent(in)    :: par, par_output
    real(real64),            intent(inout) :: grad(:)
 
    integer :: jj
    type(tdf_t) :: depsilon
    push_sub(controlfunction_gradient)
 
    select case (par%current_representation)
    case (ctr_rt)
      do jj = 1, par%dof
        ! Probably we could do without par%u, and write explicitly the values it takes.
        grad(jj) = m_two * controlfunction_alpha(par, 1) * par%u(jj, 1)*par%theta(jj) &
          - par%u(jj, 1) * tdf(par_output%f(1), jj)
      end do
    case (ctr_fourier_series, ctr_zero_fourier_series)
      do jj = 1, par%dof
        call tdf_copy(depsilon, par%f(1))
        call controlfunction_der(par, depsilon, jj)
        grad(jj) = m_two * controlfunction_alpha(par, 1) * sum(par%u(jj, :)*par%theta(:)) &
          - tdf_dot_product(par_output%f(1), depsilon)
        call tdf_end(depsilon)
      end do
 
    case (ctr_fourier_series_h, ctr_zero_fourier_series_h)
      do jj = 1, par%dof
        call tdf_copy(depsilon, par%f(1))
        call controlfunction_der(par, depsilon, jj)
        grad(jj) = - tdf_dot_product(par_output%f(1), depsilon)
        call tdf_end(depsilon)
      end do
 
    end select
 
    pop_sub(controlfunction_gradient)
  end subroutine controlfunction_gradient
  ! ---------------------------------------------------------
 
#include "controlfunction_trans_inc.F90"
 
end module controlfunction_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: