Reduced the energy grid density at the lost boundary

namelist.in

@@ -6,11 +6,11 @@
 &RUNS
  IMODE           = 1,
  NEIGEN          = 1,
- OMEGAIN(1)      = (244760.542334162, 27653.8481392209),
+ OMEGAIN(1)      = (239000, 30000),
  VSIGN_IN        = 1,
- ENERGY_IN       = 34.9,
- MUB0_IN         = 577.777777778,
- PPHI_IN         = -0.35,
+ ENERGY_IN       = 300,
+ MUB0_IN         = 300,
+ PPHI_IN         = -0.15,
  IENERGY_TPBOUND = -1,
  ENERGY_LOG_IN   = 20.,
 
@@ -25,22 +25,22 @@
  NP_COP          = 1,
  NP_CTP          = 1,
  NORBITINTSAMPLE = 2000,
- DTORBITN        = 0.5,
- DTORBITB        = 0.2,
- ERREIG          = 1.E-4,
+ DTORBITN        = 2.,
+ DTORBITB        = 1.0,
+ ERREIG          = 1.E-5,
  NMAXIT          = 25,
-&END
+&EN2
 
 ! GRID SEETINGS
 &GRID
- NRADIAL_GRID    = 30,
+ NRADIAL_GRID    = 50,
  NGTRAP_MUB0     = 10,
- TRAP_MUB0START  = 50,
+ TRAP_MUB0START  = 50.,
  TRAP_MUB0END    = 1000.,
  NGTRAP_ENERGYN  = 30,
- NGTRAP_ENERGYB  = 20,
+ NGTRAP_ENERGYB  = 15,
  TRAP_EBEND      = 20.,
- NGTRAP_PPHI     = 30 
+ NGTRAP_PPHI     = 100 
 
 &END
 
@@ -51,7 +51,7 @@
  AI = 2.0,
  R0 = 3.0,
  A  = 1.0,
- B0 = 2.8,
+ B0 = 15.,
 &END
 
 ! BULK ION PROFILES
@@ -60,7 +60,7 @@ IQ_TYPE = 1,
 ITE_TYPE = 1,
 ITI_TYPE = 1,
 TEPOLY(1) = 1.,
-TEPOLY(2) = 0.,
+TEPOLY(2) = 0.
 TEPOLY(3) = -1.,
 TIPOLY(1) = 1.,
 TIPOLY(2) = 0.,
@@ -69,16 +69,16 @@ TE0 = 3.,
 TI0 = 1.,
 QPOLY(1) = 2.,
 QPOLY(2) = 0.,
-QPOLY(3) = 2.,
+QPOLY(3) = 0.,
 &END
 
 ! FAST ION DISTRIBUTION FUNCTION
 &FAST
-NF_RATIO = 0.02
-DPPHI_NF = 0.2
-DPPHI_TPER = 0.5
-DPPHI_TPAR = 0.5
-TPER0 = 400
-TPAR0 = 100
-RRES = 2.8
+NF_RATIO = 0.005
+DPPHI_NF = 10000.
+DPPHI_TPER = 10000.
+DPPHI_TPAR = 10000.
+TPER0 = 150
+TPAR0 = 50
+RRES = 2.7
 &END

py/field.py

@@ -1,4 +1,4 @@
-from array import array
+import numpy as np
 import matplotlib.pyplot as plt
 
 f = open("field.out")
@@ -7,21 +7,23 @@
 neigen = ns[0]
 nr = ns[1]
 
-r = array('d')
-err = array('d')
-eri = array('d')
+r = np.zeros((nr), dtype=float)
+err = np.zeros((nr), dtype=float)
+eri = np.zeros((nr), dtype=float)
 
 for i1 in range(neigen):
     for i2 in range(nr):
         line = f.readline()
         number_r = map(float, line.split())
-        r.append(number_r[0])
-        err.append(number_r[1])
-        eri.append(number_r[2])
+        r[i2] = number_r[0]
+        err[i2] = (number_r[1])
+        eri[i2] = (number_r[2])
+    era = np.sqrt(np.square(err) + np.square(eri))
     plt.figure(i1)
     rep,=plt.plot(r,err)
     imp,=plt.plot(r,eri)
+    abp,=plt.plot(r,era)
     plt.ylabel('Er')
     plt.xlabel('r')
-    plt.legend([rep, imp],['Real (Mode ' + str(i1+1) + ')','Imag (Mode ' + str(i1+1) + ')'])
+    plt.legend([rep, imp, abp],['Real (Mode ' + str(i1+1) + ')','Imag (Mode ' + str(i1+1) + ')', 'Abs (Mode ' + str(i1+1) + ')'])
     plt.show()

src/paras_num.f90

@@ -32,8 +32,8 @@ module paras_num
   ! precision in finding trapped/passing boundary
   real, parameter :: errtpbound = 1e-12
   ! max iteration number in finding t/p boundary
-
-  real, parameter :: maxittpbound = 400
+  real, parameter :: maxittpbound = 1000
+
   ! maximum number of orbit harmonics
   integer, parameter :: npmax = 5
   ! maximum number of finite elements
@@ -42,7 +42,7 @@ module paras_num
   real, parameter :: ddx = 0.0001
   ! the min difference between trapedge and traplost energy to 
   ! take into account, in unit of mub0
-  real, parameter :: mindiff_el = 4e-4
+  real, parameter :: mindiff_el = 1e-4
 
   ! maximum grid points in a cubic spline interpolation
   integer, parameter :: nmaxspline = 1000

src/trap_matrix.f90

@@ -71,16 +71,22 @@ subroutine getmat3trap(this, omega, mat3)
 
           do m = this%grid(i1)%ielementmin(i2), this%grid(i1)%ielementmax(i2)
              do n = m, this%grid(i1)%ielementmax(i2)
-                do p = 1, this%grid(1)%np  
-                   tmpmat%data(n, m) = tmpmat%data(n, m) &
-                        + getint(this, i1, i2,  omega, p, m, n) &
-                        + getintnormal(this, i1, i2, -omega, p, m, n)
+                do p = 1, this%grid(1)%np
+                   if (this%grid(i1)%periodn(i2)%d(1) .eq. 0.) then
+                      ! fqc constant for given mub0 and pphi, get int will fail
+                      tmpmat%data(n, m) = tmpmat%data(n, m) &
+                           + getintnormal(this, i1, i2,  omega, p, m, n) &
+                           + getintnormal(this, i1, i2, -omega, p, m, n)
+                   else
+                      tmpmat%data(n, m) = tmpmat%data(n, m) &
+                           + getint(this, i1, i2,  omega, p, m, n) &
+                           + getintnormal(this, i1, i2, -omega, p, m, n)
+                   end if
                 end do
                 tmpmat%data(m, n) = tmpmat%data(n, m)
              end do
           end do
 
-
           if (i2 .eq. 1) then
              dpphi = this%grid(i1)%pphigrid(2) - this%grid(i1)%pphigrid(1)
              tmpmat2 = cmplx(-0.5*dpphi) * tmpmat
@@ -240,6 +246,7 @@ complex function getint(this, imub0, ipphi, omega, p, m, n)
     ! OUTPUT: the value of the integral
 
     use cubic, only : getcoeffs
+    use distribution_fun
     implicit none
 
     type(tmatrix), intent(in) :: this
@@ -251,6 +258,7 @@ complex function getint(this, imub0, ipphi, omega, p, m, n)
     integer :: i1, i2, iup, ipphib
     logical :: ltype1
     real, dimension(4) :: c, pc
+    real :: mub0, pphi, ee
 
     periodp = 2 * pi * real(p) / omega
 
@@ -303,6 +311,7 @@ complex function getint(this, imub0, ipphi, omega, p, m, n)
        pc(3) = this%grid(imub0)%periodn(ipphi)%b(i1)
        pc(4) = this%grid(imub0)%periodn(ipphi)%a(i1)
        tmpres = landauint(c, pc, periodp, x1, x2)
+
        results = results + tmpres
     end do
     getint = results