text:
      SUBROUTINE GETVEC( IPTR, X, NMAX, N )
      INTEGER*4       IPTR, NMAX, N
      REAL*8          X( NMAX )
C
C  Copy input argument pointed to by IPTR into the target vector X.
C  Make sure the input argument is a real*8 vector, and that
C  it doesn't exceed the length of the target vector.
C  Put the length of the input vector in N.
C
C  MEX functions
C
      INTEGER*4       MXGETPR
      INTEGER*4       MXGETM, MXGETN
      INTEGER*4       MXISDOUBLE
C
C  Local variables
C
      INTEGER*4       NROWS, NCOLS, M, PR
      CHARACTER*100   MSG
C
      NROWS = MXGETM( IPTR )
      NCOLS = MXGETN( IPTR )
      M     = MIN( NROWS, NCOLS )
      N     = MAX( NROWS, NCOLS )
      IF ( M .NE. 1 ) THEN
         MSG = 'INPUT ARGUMENT MUST BE A VECTOR.'
         CALL MEXERRMSGTXT( MSG )
      END IF
      IF ( N .GT. NMAX ) THEN
         MSG = 'INPUT ARGUMENT IS TOO LONG.'
         CALL MEXERRMSGTXT( MSG )
      END IF
      IF ( MXISDOUBLE( IPTR ) .EQ. 0 ) THEN
         MSG = 'INPUT ARGUMENT MUST BE REAL*8.'
         CALL MEXERRMSGTXT( MSG )
      END IF
      PR = MXGETPR( IPTR )
      CALL MXCOPYPTRTOREAL8( PR, X, N )
      RETURN
      END
C
      SUBROUTINE CNVSPR( NCOL, NNZSH, SH, IRNSH, ICNSH, A, IROW, JCOL )
      INTEGER*4       NCOL, NNZSH
      INTEGER*4       IRNSH( NNZSH ), ICNSH( NNZSH ), IROW( NNZSH )
      INTEGER*4       JCOL( NCOL + 1 )
      REAL*8          SH( NNZSH ), A( NNZSH )
C
C  Convert the CUTE sparse matrix described by SH, IRNSH and ICNSH
C  to the MATLAB sparse matrix described by A, IROW and JCOL.
C
C  Local variables
C
      INTEGER*4       NCOLP1, I, J, K, JJ
C
C  Zero out JCOL and then use it to store the number of nonzero
C  entries in column J of the sparse matrix.
C
      NCOLP1 = NCOL + 1
      DO 100 J = 1, NCOLP1
         JCOL( J ) = 0
  100 CONTINUE
      DO 110 K = 1, NNZSH
         J         = ICNSH( K )
         JCOL( J ) = JCOL( J ) + 1
  110 CONTINUE
      JCOL( NCOLP1 ) = NNZSH + 1
C
C  Now go backwards through JCOL to find the starting index for the
C  nonzero entries in column J.
C
      DO 120 J = NCOL, 1, -1
         JCOL( J ) = JCOL( J + 1 ) - JCOL( J )
  120 CONTINUE
C
C  Copy the entries from SH and IRNSH into the correct entries
C  in A and IROW.  Since MATLAB is written in C, subtract 1 from
C  each of the entries in IROW.  Use JCOL to keep track of what
C  place we're in for each column.
C
      DO 130 K = 1, NNZSH
         I          = IRNSH( K )
         J          = ICNSH( K )
         JJ         = JCOL ( J )
         A   ( JJ ) = SH   ( K )
         IROW( JJ ) = I  - 1
         JCOL( J  ) = JJ + 1
  130 CONTINUE
C
C  Restore the entries in JCOL to the starting index for the nonzero
C  entries in each column.  Since MATLAB is written in C, subtract 1
C  from each of the entries in JCOL.
C
      DO 140 J = NCOL, 2, -1
         JCOL( J ) = JCOL( J - 1 ) - 1
  140 CONTINUE
      JCOL( 1 ) = 0
      JCOL( NCOLP1 ) = JCOL( NCOLP1 ) - 1
      RETURN
      END
C
      SUBROUTINE CNVSAB( NCOL, NNZSH, SH, IRNSH, ICNSH, NNZA, A, IROWA,
     *                   JCOLA, NNZB, B, IROWB, JCOLB )
      INTEGER*4       NCOL, NNZSH, NNZA, NNZB
      INTEGER*4       IRNSH( NNZSH ), ICNSH( NNZSH ), IROWA( NNZSH )
      INTEGER*4       IROWB( NCOL )
      INTEGER*4       JCOLA( NCOL + 1 ), JCOLB( NCOL + 1 )
      REAL*8          SH( NNZSH ), A( NNZSH ), B( NCOL )
C
C  The CUTE sparse matrix described by SH, IRNSH and ICNSH contains
C  both a sparse matrix and a sparse vector.  The k-th nonzero entry
C  in SH belongs to the sparse vector if IRNSH(k)=0, and belongs
C  to row i of the sparse matrix if IRNSH(k)=i>0.
C  This subroutine returns the MATLAB sparse matrix described by
C  NNZA, A, IROWA, JCOLA, and the MATLAB sparse vector described by
C  NNZB, B, IROWB, JCOLB.
C
C  Local variables
C
      INTEGER*4       NCOLP1, I, J, K, JJ
C
C  Zero out JCOLA and JCOLB.  Use JCOLA to store the number of nonzero
C  entries in column J of the sparse matrix, and use JCOLB to store
C  the number of nonzero entries in column J of the sparse vector.
C
      NCOLP1 = NCOL + 1
      DO 100 J = 1, NCOLP1
         JCOLA( J ) = 0
         JCOLB( J ) = 0
  100 CONTINUE
      NNZA = 0
      NNZB = 0
      DO 110 K = 1, NNZSH
         I         = IRNSH( K )
         J         = ICNSH( K )
         IF ( I .GT. 0 ) THEN
            JCOLA( J ) = JCOLA( J ) + 1
            NNZA       = NNZA + 1
         ELSE
            JCOLB( J ) = JCOLB( J ) + 1
            NNZB       = NNZB + 1
         END IF
  110 CONTINUE
      JCOLA( NCOLP1 ) = NNZA + 1
      JCOLB( NCOLP1 ) = NNZB + 1
C
C  Now go backwards through JCOLA and JCOLB to find the starting index
C  for the nonzero entries in column J.
C
      DO 120 J = NCOL, 1, -1
         JCOLA( J ) = JCOLA( J + 1 ) - JCOLA( J )
         JCOLB( J ) = JCOLB( J + 1 ) - JCOLB( J )
  120 CONTINUE
C
C  Copy the entries from SH and IRNSH into the correct entries
C  in A and IROWA, or B if IRNSH(k)=0.  Since MATLAB is written
C  in C, subtract 1 from each of the entries in IROWA.  Set all
C  of the entries in IROWB=0.  Use JCOLA and JCOLB to keep track
C  of what place we're in for each column.
C
      DO 130 K = 1, NNZSH
         I          = IRNSH( K )
         J          = ICNSH( K )
         IF ( I .GT. 0 ) THEN
            JJ          = JCOLA( J )
            A    ( JJ ) = SH   ( K )
            IROWA( JJ ) = I  - 1
            JCOLA( J  ) = JJ + 1
         ELSE
            JJ          = JCOLB( J )
            B    ( JJ ) = SH   ( K )
            IROWB( JJ ) = 0
            JCOLB( J  ) = JJ + 1
         END IF
  130 CONTINUE
C
C  Restore the entries in JCOLA and JCOLB to the starting index for
C  the nonzero entries in each column.  Since MATLAB is written in C,
C  subtract 1 from each of the entries in JCOLA and JCOLB.
C
      DO 140 J = NCOL, 2, -1
         JCOLA( J ) = JCOLA( J - 1 ) - 1
         JCOLB( J ) = JCOLB( J - 1 ) - 1
  140 CONTINUE
      JCOLA( 1 ) = 0
      JCOLB( 1 ) = 0
      JCOLA( NCOLP1 ) = JCOLA( NCOLP1 ) - 1
      JCOLB( NCOLP1 ) = JCOLB( NCOLP1 ) - 1
      RETURN
      END
C
      SUBROUTINE CNVCHR( CNAMES, NWORD, LWORD, WORK )
      CHARACTER*(*) CNAMES( NWORD )
      INTEGER*4     NWORD, LWORD
      REAL*8        WORK( NWORD, LWORD )
C
C  Convert individual character entries in character array CNAMES
C  to the floating point values of their position in the collating
C  sequence (i.e., convert to floating point ASCII values).
C  Put these values in the REAL*8 array WORK.
C
      DO 120 I = 1, NWORD
         DO 100 J = 1, LWORD
            WORK( I, J ) = DFLOAT( ICHAR( CNAMES(I)(J:J) ) )
  100    CONTINUE
  120 CONTINUE
      RETURN
      END
C
C     UPPER - CONVERT STRING TO UPPER CASE EXCEPT WHERE ENCLOSED IN QUOTES
C
C     SOURCE:
C     FORTRAN TOOLS FOR VAX/VMS AND MS-DOS
C     BY RUSSELL K. JONES AND TRACY CRABTREE,
C     JOHN WILEY & SONS, 1988.
C
      CHARACTER*(*) FUNCTION UPPER(STRING)
      CHARACTER*(*) STRING
      CHARACTER*1   NULL, QUOTE
      INTEGER       I
C
      NULL = CHAR(0)
      QUOTE = CHAR(39)
      I = 1
      DO WHILE (STRING(I:I) .NE. NULL)
         IF (STRING(I:I) .EQ. QUOTE) THEN
            UPPER(I:I) = QUOTE
            I = I + 1
            DO WHILE (STRING(I:I) .NE. QUOTE .AND. STRING(I:I) .NE.
     &       NULL)
               UPPER(I:I) = STRING(I:I)
               I = I + 1
            END DO
            UPPER(I:I) = STRING(I:I)
            IF (STRING(I:I) .NE. NULL) I = I + 1
         ELSE
            IF (LLT(STRING(I:I),'a') .OR. LGT(STRING(I:I),'z')) THEN
               UPPER(I:I) = STRING(I:I)
            ELSE
               UPPER(I:I) =
     &          CHAR(ICHAR('A') - ICHAR('a') + ICHAR(STRING(I:I)))
            END IF
            I = I + 1
         END IF
      END DO
      UPPER(I:I) = NULL
      RETURN
      END
C
C     NULLSTR - NULL-TERMINATE AN UNTERMINATED STRING
C
C     SOURCE:
C     FORTRAN TOOLS FOR VAX/VMS AND MS-DOS
C     BY RUSSELL K. JONES AND TRACY CRABTREE,
C     JOHN WILEY & SONS, 1988.
C
      CHARACTER*(*) FUNCTION NULLSTR(STR)
      CHARACTER*(*) STR
      CHARACTER*1   NULL, BLANK
      INTEGER       I
      NULL = CHAR(0)
      BLANK = CHAR(32)
      I = LEN(STR)
      DO WHILE (I .GT. 0 .AND. STR(I:I) .EQ. BLANK)
         I = I - 1
      END DO
      IF (I .EQ. 0) THEN
         NULLSTR = NULL
      ELSE
         NULLSTR = STR(1:I)//NULL
      END IF
      RETURN
      END
C
C     EQUAL - COMPARE TWO STRINGS FOR EQUALITY
C
C     SOURCE:
C     FORTRAN TOOLS FOR VAX/VMS AND MS-DOS
C     BY RUSSELL K. JONES AND TRACY CRABTREE,
C     JOHN WILEY & SONS, 1988.
C
      LOGICAL FUNCTION EQUAL(S,T)
      CHARACTER*(*) S, T
      CHARACTER*1   NULL
      INTEGER       I
      NULL = CHAR(0)
      I = 1
      DO WHILE (S(I:I) .NE. NULL)
         IF (S(I:I) .NE. T(I:I)) THEN
            EQUAL = .FALSE.
            RETURN
         END IF
         I = I + 1
      END DO
      IF (T(I:I) .EQ. NULL) THEN
         EQUAL = .TRUE.
      ELSE
         EQUAL = .FALSE.
      END IF
      RETURN
      END
C

max_forks_repo_name: dynamics-of-stellar-systems/dynamite


text:
! Written by Ulf Saalmann
MODULE spaic2_betaSpect
IMPLICIT NONE

  ! GLOBAL PARAMETERS
  include "spaic2_defs.f90"

  integer :: num_qns

  double precision :: freqmin=0.005d0, freqmax=0.080d0
  integer, parameter :: nfreq=75
  double precision :: eb

  ! ALLOCATABLES
  double precision, allocatable :: pots(:,:)
  double precision, allocatable :: rr(:)
  double precision, allocatable :: bpot(:)
  double precision, allocatable :: cpot(:)
  double precision, allocatable :: bpsi(:)
  double precision, allocatable :: cpsi(:)
  double precision, allocatable :: pot(:)
  double precision, allocatable :: pot0(:)
  double precision, allocatable :: pot1(:)
  double precision, allocatable :: spec(:, :)
  double precision, allocatable :: beta(:, :)
  double precision, allocatable :: potcha(:)
  double precision, allocatable :: rhos(:,:)
  double precision, allocatable :: freq(:)
  double precision, allocatable :: para(:)
  integer, allocatable :: quantum_numbers(:, :)

  ! INTERFACE

  ! SUBROUTINES
  public :: initialize
  public :: crossection
  public :: deinitialize
  public :: compute_beta
  public :: free_all_memory
  public :: set_density_and_qns

CONTAINS

  SUBROUTINE free_all_memory()
  IMPLICIT NONE
    call deinitialize()
    if(allocated(para)) then
      deallocate(para)
    endif
    if(allocated(quantum_numbers)) then
      deallocate(quantum_numbers)
    endif
  END SUBROUTINE free_all_memory

  SUBROUTINE set_density_and_qns(params, params_size, qns, qns_size)
    integer,                                  intent(in) :: params_size
    double precision, dimension(params_size), intent(in) :: params
    integer,                                  intent(in) :: qns_size
    integer,          dimension(qns_size, 2), intent(in) :: qns

    ! para = params
    if(allocated(para)) then
      deallocate(para)
    endif
    allocate(para(1:params_size))
    para(:) = params(:)
    num_para = params_size

    ! quantum_numbers = qns
    if(allocated(quantum_numbers)) then
      deallocate(quantum_numbers)
    endif
    allocate(quantum_numbers(1:qns_size, 1:2))
    quantum_numbers(:, :) = qns(:, :)
    num_qns = qns_size

    call initialize()
  END SUBROUTINE set_density_and_qns

  SUBROUTINE compute_beta()
    integer :: ii, l, n
    double precision :: econt(0:nfreq), dipm(0:nfreq), dipp(0:nfreq), &
      phim(0:nfreq), phip(0:nfreq)

    ! calculate differential cross section, beta parameter, etc.
    do ii=1,num_qns
      l = quantum_numbers(ii, 1)
      n = quantum_numbers(ii, 2)
      if (l==0) then
        dipm(:) = 0
        call crossection(l, n, l+1, freq, nfreq, econt,dipp,phip)
      else
        call crossection(l, n, l-1, freq, nfreq, econt,dipm,phim)
        call crossection(l, n, l+1, freq, nfreq, econt,dipp,phip)
      end if

      beta(:, ii) = (l*(l-1) * dipm(:)**2 + (l+1)*(l+2) * dipp(:)**2 &
           - 6*l*(l+1) * dipm(:)*dipp(:) * cos(phip(:)-phim(:))) &
           / ((2*l+1) * (l*dipm(:)**2 + (l+1)*dipp(:)**2) + 1d-99)
      spec(:, ii) = dipm(:)**2 + dipp(:)**2

    end do

  END SUBROUTINE compute_beta

  include "spaic2_misc.f90"
  include "spaic2_cwf.f90"

  SUBROUTINE crossection(lin, nin, lfin, freq, nfreq, econt, dip, phi)
  IMPLICIT NONE
    integer, intent(in) :: lin,nin,lfin, nfreq
    double precision, intent(in) :: freq(0:nfreq)
    double precision, intent(out) :: econt(0:nfreq),dip(0:nfreq),phi(0:nfreq)
    integer :: i, j
    double precision :: kc, angular,s,c

    j = max(lin,lfin)
    angular = j/sqrt(4*j*j-1d0)

    bpot(:) = pot(:) + lin*(lin+1)*pot1(:)
    cpot(:) = pot(:) + lfin*(lfin+1)*pot1(:)

    ! calculate bound state
    call bdiag(nin, bpot(:), nr,eb,bpsi)

    ! calculate continum states and dipole matrix element
    do i=0, nfreq
      econt(i) = freq(i)
      kc = sqrt(2*econt(i))
      if (econt(i)<ecmin) then
        dip(i) = 0
        phi(i) = 0
        cpsi(:) = 0
        s = 1
        c = 0
      else
        call cnumerov(lfin,econt(i),cpot(0:), nr,cpsi,s,c)
        !! write(0,'(3i1,i4,f8.3,2f8.4)') lin,nin,lfin,i,27.21*econt(i),s,c
        phi(i) = atan2(c,s)
        dip(i) = angular * sum(bpsi(:)*rr(:)*cpsi(:))/sqrt(s**2+c**2)
      end if

    end do

  END SUBROUTINE crossection

  subroutine initialize()
    implicit none
    integer :: i,ii
    double precision :: r,drbspline

    call deinitialize()

    allocate(pots(0:nr, 1:num_para))
    allocate(rr(0:nr), bpot(0:nr), cpot(0:nr), bpsi(0:nr), cpsi(0:nr))
    allocate(pot(0:nr), pot0(0:nr), pot1(0:nr))
    allocate(spec(0:nfreq, 1:num_qns), beta(0:nfreq, 1:num_qns))
    allocate(potcha(1:num_para), rhos(0:nr, 1:num_para))
    allocate(freq(0:nfreq))

    ! Compute frequency axis
    do i = 0, nfreq
      freq(i) = freqmin + i * ((freqmax-freqmin)/nfreq)
    end do

    ! calculate all potentials
    do ii=0,nr
      r = ii*dr + 1d-6
      rr(ii) = r
      if (ii>0) pot1(ii) = 1/(2*r**2)
    end do
    pot1(0) = pot1(1)

    ! make the knot distance that the last b-spline vanishes at rpot
    drbspline = rpot/(num_para+2)
    do ii = 0, nr
      r = rr(ii) + 1d-10
      pot0(ii) =0
      do i=1,num_para
        pots(ii,i) = bspline3intnorm(i,r,drbspline)
      end do
    end do

    pot0(:) = 0
    pot(:)  = 0
    do i = 1, num_para
      ! shift the parameter's c.o.m. to guarantee sum(para(:)) = 1
      pot(:) = pot(:) + (para(i)-(sum(para(:))-1)/num_para)*pots(:,i)
      ! scale all parameters to guarantee sum(para(:)) = 1
      ! pot(:) = pot(:) + (para(i)/sum(para(:)))*pots(:,i)
    end do

  END SUBROUTINE initialize

  SUBROUTINE deinitialize()
  IMPLICIT NONE
    if(allocated(pots)) then
      deallocate(pots)
    endif
    if(allocated(rr)) then
      deallocate(rr)
    endif
    if(allocated(bpot)) then
      deallocate(bpot)
    endif
    if(allocated(cpot)) then
      deallocate(cpot)
    endif
    if(allocated(bpsi)) then
      deallocate(bpsi)
    endif
    if(allocated(cpsi)) then
      deallocate(cpsi)
    endif
    if(allocated(pot)) then
      deallocate(pot)
    endif
    if(allocated(pot0)) then
      deallocate(pot0)
    endif
    if(allocated(pot1)) then
      deallocate(pot1)
    endif
    if(allocated(spec)) then
      deallocate(spec)
    endif
    if(allocated(beta)) then
      deallocate(beta)
    endif
    if(allocated(potcha)) then
      deallocate(potcha)
    endif
    if(allocated(rhos)) then
      deallocate(rhos)
    endif
    if(allocated(freq)) then
      deallocate(freq)
    endif
  END SUBROUTINE deinitialize

END MODULE spaic2_betaSpect

max_stars_repo_name: iaaft/spyclib

max_issues_repo_name: iaaft/spyclib

max_forks_repo_name: iaaft/spyclib