istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

add a switch for realtime mode

Browse Source 
This adds a single switch to configure GRAY for realtime operation.

In realtime a single ray is traced until absorption reaches ≈50% of the
total. The simulation is stopped immediately after returning only the
position (as ρ_p=√ψ) of the peak. So, no post-processing is performed at
all. In addition:

 - all outputs units are inactivated (equivalent to passing --units 0);

 - current drive computation is disabled (ecrh_cd.ieccd=0);

 - absorption is computed in weakly relativistic approximation
   (ecrh_cd.iwarm=1);

 - reflections and multiple plasma passes are disabled;
This commit is contained in:
Michele Guerini Rocco 2022-05-24 11:19:30 +02:00
parent 92b3ad9bd3
commit 7342566ac0
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
7 changed files with 218 additions and 76 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
input/gray.ini
View File

@ -39,6 +39,18 @@ integrator = 4
; the initial step size.
adaptive_step = false
; Enable the realtime mode
; In realtime a single ray is traced until absorption reaches ≈50% of the
; total. The simulation is stopped immediately after returning only the
; position (as ρ_p=√ψ) of the peak. So, no post-processing is
; performed at all. In addition:
; - all outputs units are inactivated (equivalent to passing --units 0);
; - current drive computation is disabled (ecrh_cd.ieccd=0);
; - absorption is computed in weakly relativistic approximation
; (ecrh_cd.iwarm=1);
; - reflections and multiple plasma passes are disabled;
realtime = false
[ecrh_cd]
; Choice of the power absorption model

21
src/beamdata.f90
View File

@ -8,7 +8,7 @@ module beamdata
contains
subroutine init_btr(rtrparam,ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
tau0,alphaabs0,dpds,dids0,ccci0,p0jk,ext,eyt,iiv)
use gray_params, only : raytracing_parameters
use const_and_precisions, only : two
implicit none
@ -17,7 +17,7 @@ contains
gri,psjki,ppabs,ccci
real(wp_), dimension(:,:,:), intent(out), pointer :: xc,du1,ggri
real(wp_), dimension(:), intent(out), pointer :: tau0,alphaabs0, &
dids0,ccci0,p0jk
dpds,dids0,ccci0,p0jk
complex(wp_), dimension(:), intent(out), pointer :: ext, eyt
integer, dimension(:), intent(out), pointer :: iiv
@ -52,7 +52,7 @@ contains
! yp = dy/dt,
! etc.
call alloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
tau0,alphaabs0,dpds,dids0,ccci0,p0jk,ext,eyt,iiv)
end subroutine init_btr
@ -166,36 +166,36 @@ contains
subroutine alloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
tau0,alphaabs0,dpds,dids0,ccci0,p0jk,ext,eyt,iiv)
implicit none
real(wp_), dimension(:,:), intent(out), pointer :: ywork,ypwork, &
gri,psjki,ppabs,ccci
real(wp_), dimension(:,:,:), intent(out), pointer :: xc,du1,ggri
real(wp_), dimension(:), intent(out), pointer :: tau0,alphaabs0, &
dids0,ccci0,p0jk
dpds,dids0,ccci0,p0jk
complex(wp_), dimension(:), intent(out), pointer :: ext, eyt
integer, dimension(:), intent(out), pointer :: iiv
call dealloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
tau0,alphaabs0,dpds,dids0,ccci0,p0jk,ext,eyt,iiv)
allocate(ywork(6,nray), ypwork(6,nray), gri(3,nray), ggri(3,3,nray), &
xc(3,nrayth,nrayr), du1(3,nrayth,nrayr), &
psjki(nray,nstep), ppabs(nray,nstep), ccci(nray,nstep), &
tau0(nray), alphaabs0(nray), dids0(nray), ccci0(nray), &
p0jk(nray), ext(nray), eyt(nray), iiv(nray))
tau0(nray), alphaabs0(nray), dpds(nstep), dids0(nray), &
ccci0(nray), p0jk(nray), ext(nray), eyt(nray), iiv(nray))
end subroutine alloc_beam
subroutine dealloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
tau0,alphaabs0,dpds,dids0,ccci0,p0jk,ext,eyt,iiv)
implicit none
real(wp_), dimension(:,:), intent(out), pointer :: ywork,ypwork, &
gri,psjki,ppabs,ccci
real(wp_), dimension(:,:,:), intent(out), pointer :: xc,du1,ggri
real(wp_), dimension(:), intent(out), pointer :: tau0,alphaabs0, &
dids0,ccci0,p0jk
dpds,dids0,ccci0,p0jk
complex(wp_), dimension(:), intent(out), pointer :: ext, eyt
integer, dimension(:), intent(out), pointer :: iiv
@ -210,6 +210,7 @@ contains
if (associated(ccci)) deallocate(ccci)
if (associated(tau0)) deallocate(tau0)
if (associated(alphaabs0)) deallocate(alphaabs0)
if (associated(dpds)) deallocate(dpds)
if (associated(dids0)) deallocate(dids0)
if (associated(ccci0)) deallocate(ccci0)
if (associated(p0jk)) deallocate(p0jk)

121
src/gray_core.f90
View File

@ -82,8 +82,10 @@ contains
real(wp_), dimension(:,:), pointer :: psjki=>null(),ppabs=>null(),ccci=>null()
real(wp_), dimension(:,:), pointer :: taus=>null(),stnext=>null(), &
yw0=>null(),ypw0=>null(),cpls=>null()
! Note: dpds is only used in realtime mode to compute the power absorption peak
real(wp_), dimension(:), pointer :: p0ray=>null(),tau0=>null(),alphaabs0=>null(), &
dids0=>null(),ccci0=>null(),tau1=>null(),etau1=>null(),cpl1=>null(),lgcpl1=>null()
dpds=>null(), dids0=>null(),ccci0=>null(),tau1=>null(),etau1=>null(),cpl1=>null(), &
lgcpl1=>null()
real(wp_), dimension(:), pointer :: p0jk=>null()
real(wp_), dimension(:), pointer :: jphi_beam=>null(),pins_beam=>null(), &
currins_beam=>null(), dpdv_beam=>null(),jcd_beam=>null(),stv=>null(), &
@ -109,32 +111,46 @@ contains
! Initialise the ray variables (beamtracing)
call init_btr(params%raytracing, yw, ypw, xc, du1, gri, ggri, psjki, ppabs, ccci, &
tau0, alphaabs0, dids0, ccci0, p0jk, ext, eyt, iiv)
tau0, alphaabs0, dpds, dids0, ccci0, p0jk, ext, eyt, iiv)
! Initialise the dispersion module
if (params%ecrh_cd%iwarm > 1) call expinit
! In realtime mode the raytracing is stopped as soon as
! possible and all these extra features are unnecessary
if (params%raytracing%realtime) then
nnd = 0
else
! Initialise the magsurf_data module
call flux_average ! requires frhotor for dadrhot,dvdrhot
! Note: 1. needed for computing flux surface averages (dP/dV, etc.)
! 2. requires frhotor for dadrhot, dvdrhot
call flux_average
! Initialise the output profiles
! Initialise output profiles module
call pec_init(params%output%ipec, rhout)
nnd = size(rhop_tab) ! number of radial profile points
! Allocate memory for the results...
allocate(results%dpdv(params%output%nrho))
allocate(results%jcd(params%output%nrho))
! ...and initialise them
results%pabs = zero
results%icd = zero
results%dpdv = zero
results%jcd = zero
end if
! Initialise multipass module
call alloc_multipass(nnd, iwait, iroff, iop, iow, yynext, yypnext, &
yw0, ypw0, stnext, stv, dst, p0ray, taus, tau1, &
etau1, cpls, cpl1, lgcpl1, jphi_beam, pins_beam, &
currins_beam, dpdv_beam, jcd_beam, psipv, chipv)
! ...and initialise them
results%pabs = zero
results%icd = zero
results%dpdv = zero
results%jcd = zero
! ========= set environment END =========
! ======== pre-proc prints BEGIN ========
if (.not. params%raytracing%realtime) then
call print_headers(params)
! print ψ surface for q=1.5 and q=2 on file and log psi,rhot,rhop
@ -155,6 +171,7 @@ contains
call print_maps(bres, xgcn, &
norm2(params%antenna%pos(1:2)) * 0.01_wp_, &
sin(params%antenna%beta*degree))
end if
! ========= pre-proc prints END =========
! =========== main loop BEGIN ===========
@ -220,7 +237,7 @@ contains
cycle
end if
call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dpds,dids0,ccci0,iiv)
if(ip == 1) then ! 1st pass
igrad_b = params%raytracing%igrad ! * input value, igrad_b=0 from 2nd pass
@ -240,8 +257,11 @@ contains
zzm = yw(3,jk)*0.01_wp_
rrm = sqrt(yw(1,jk)*yw(1,jk)+yw(2,jk)*yw(2,jk))*0.01_wp_
if(inside(data%equilibrium%rlim, data%equilibrium%zlim, &
nlim, rrm, zzm)) then ! * start propagation in/outside vessel?
! Note: `inside` call is expensive and can be skipped in realtime
if (params%raytracing%realtime) then
iow(jk) = 1
else if (inside(data%equilibrium%rlim, data%equilibrium%zlim, & ! * start propagation in/outside vessel?
nlim, rrm, zzm)) then
iow(jk) = 1 ! + inside
else
iow(jk) = 0 ! + outside
@ -305,9 +325,13 @@ contains
zzm = xv(3)*0.01_wp_
rrm = sqrt(xv(1)*xv(1)+xv(2)*xv(2))*0.01_wp_
! Note: `inside` call is expensive and can be skipped in realtime
if (params%raytracing%realtime) then
ins_wl = .true.
else
ins_wl = inside(data%equilibrium%rlim, data%equilibrium%zlim, nlim,rrm,zzm) ! in/out vessel?
end if
ins_pl = (psinv>=zero .and. psinv<params%profiles%psnbnd) ! in/out plasma?
ins_wl = inside(data%equilibrium%rlim, data%equilibrium%zlim, &
nlim,rrm,zzm) ! in/out vessel?
ent_pl = (mod(iop(jk),2) == 0 .and. ins_pl) ! enter plasma
ext_pl = (mod(iop(jk),2) == 1 .and. .not.ins_pl) ! exit plasma
ent_wl = (mod(iow(jk),2) == 0 .and. ins_wl) ! enter vessel
@ -474,7 +498,8 @@ contains
pow = p0ray(jk) * exp(-tau) ! residual power: P = P₀exp(-τ)
ppabs(jk,i) = p0ray(jk) - pow ! absorbed power: P_abs = P₀ - P
dids = didp * pow * alpha ! current driven: dI/ds = dI/dP⋅dP/ds = dI/dP⋅P⋅α
dpds(i) = pow * alpha ! power density: dP/ds = P⋅α
dids = didp * dpds(i) ! current driven: dI/ds = dI/dP⋅dP/ds
! current: I = ∫dI/ds⋅ds using the trapezoidal rule
ccci(jk,i) = ccci0(jk) + 0.5_wp_*(dids0(jk) + dids) * dersdst * dst(jk)
@ -503,6 +528,26 @@ contains
end do
end if
if (params%raytracing%realtime .and. pow <= 0.4_wp_*p0ray(1)) then
! Compute the approximate position of the absorption peak
block
use utils, only : vmax, parabola_vertex
real(wp_) :: power, peak(2)
integer :: index
! Find maximum in dP/ds
call vmax(dpds, i, power, index)
! Interpolate the peak with a parabola
peak = parabola_vertex(psjki(1, index - 1:index + 1), &
dpds(index - 1:index + 1))
results%rho_peak = sqrt(peak(1))
end block
! Stop propagation loop
exit
end if
! print ray positions for j=nrayr in local reference system
if(mod(i,params%output%istpr) == 0) then
if(params%raytracing%nray > 1 .and. all(.not.iwait)) &
@ -533,9 +578,11 @@ contains
icd_beam = sum(ccci(:,i))
call vmaxmin(tau0,params%raytracing%nray,taumn,taumx) ! taumn,taumx for print
! compute power and current density profiles for all rays
call spec(psjki,ppabs,ccci,iiv,pabs_beam,icd_beam,dpdv_beam,jphi_beam,jcd_beam, &
pins_beam,currins_beam)
if (.not. params%raytracing%realtime) then
! Compute power and current density profiles for all rays
call spec(psjki, ppabs, ccci, iiv, pabs_beam, icd_beam, &
dpdv_beam, jphi_beam, jcd_beam, pins_beam, currins_beam)
end if
pabs_pass(iox) = pabs_pass(iox) + pabs_beam ! 0D results for current pass, sum on O/X mode beams
icd_pass(iox) = icd_pass(iox) + icd_beam
@ -578,19 +625,32 @@ contains
end if
end if
call print_pec(rhop_tab,rhot_tab,jphi_beam,jcd_beam,dpdv_beam,currins_beam, &
pins_beam,ip) ! *print power and current density profiles for current beam
if (.not. params%raytracing%realtime) then
call print_pec(rhop_tab, rhot_tab, jphi_beam, jcd_beam, &
dpdv_beam, currins_beam, pins_beam, ip) ! *print power and current density profiles for current beam
call postproc_profiles(pabs_beam,icd_beam,rhot_tab,dpdv_beam,jphi_beam, &
rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip,rhotp,drhotp,rhotj, &
drhotj,dpdvmx,jphimx,ratjamx,ratjbmx) ! *compute profiles width for current beam
call postproc_profiles(pabs_beam, icd_beam, rhot_tab, dpdv_beam, &
jphi_beam, rhotpav, drhotpav, rhotjava, &
drhotjava, dpdvp, jphip, rhotp, drhotp, &
rhotj, drhotj, dpdvmx, jphimx, ratjamx, ratjbmx) ! *compute profiles width for current beam
call print_finals(pabs_beam,icd_beam,dpdvp,jphip,rhotpav,rhotjava, &
drhotpav,drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx, &
ratjbmx,stv(1),psipv(index_rt),chipv(index_rt),index_rt,sum(p0ray), &
cpl_beam1,cpl_beam2) ! *print 0D results for current beam
call print_finals(pabs_beam, icd_beam, dpdvp, jphip, rhotpav, &
rhotjava, drhotpav, drhotjava, dpdvmx, jphimx, &
rhotp, rhotj, drhotp, drhotj, ratjamx, ratjbmx, &
stv(1), psipv(index_rt), chipv(index_rt), &
index_rt, sum(p0ray), cpl_beam1, cpl_beam2) ! *print 0D results for current beam
end if
! ============ post-proc END ============
! Store the accurate position of the absorption peak
if (.not. params%raytracing%realtime .and. ip == 1 .and. rhotpav /= 0) then
! Note: rho_peak is poloidal, we covert ρ_t to ρ_p
block
use equilibrium, only: frhopol
results%rho_peak = frhopol(rhotpav)
end block
end if
end do beam_loop
call log_debug('beam loop end', mod='gray_core', proc='gray_main')
! ============ beam loop END ============
@ -619,7 +679,7 @@ contains
! ========== free memory BEGIN ==========
call dealloc_surfvec
call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci,tau0, &
alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
alphaabs0,dpds,dids0,ccci0,p0jk,ext,eyt,iiv)
call dealloc_pec
call dealloc_multipass(iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0, &
stnext,stv,dst,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, &
@ -671,12 +731,12 @@ contains
subroutine vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
subroutine vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dpds,dids0,ccci0,iiv)
use const_and_precisions, only : zero
implicit none
! arguments
real(wp_), dimension(:,:), intent(out) :: psjki,ppabs,ccci
real(wp_), dimension(:), intent(out) :: tau0,alphaabs0,dids0,ccci0
real(wp_), dimension(:), intent(out) :: tau0,alphaabs0,dpds,dids0,ccci0
integer, dimension(:), intent(out) :: iiv
!! common/external functions/variables
! integer :: jclosest
@ -693,6 +753,7 @@ contains
ccci = zero
tau0 = zero
alphaabs0 = zero
dpds = zero
dids0 = zero
ccci0 = zero
iiv = 1

6
src/gray_params.f90
View File

@ -62,6 +62,7 @@ module gray_params
logical :: adaptive_step ! Allow variable step sizes
integer :: ipass ! Number of plasma passes
integer :: ipol ! Whether to compute wave polarisation
logical :: realtime ! Enable the realtime mode
end type
! EC resonant heating & Current Drive parameters
@ -132,6 +133,7 @@ module gray_params
type gray_results
real(wp_) :: pabs ! Total absorbed power
real(wp_) :: icd ! Total driven current
real(wp_) :: rho_peak ! Position of the absoprtion peak
real(wp_), allocatable :: dpdv(:) ! Absorbed power density
real(wp_), allocatable :: jcd(:) ! Driven current density
end type
@ -375,6 +377,7 @@ contains
! gray_params.data has been deprecated
params%raytracing%integrator = 4
params%raytracing%adaptive_step = .false.
params%raytracing%realtime = .false.
end subroutine read_gray_params
@ -397,7 +400,8 @@ contains
istpr0 = params%output%istpr
istpl0 = params%output%istpl
if (params%raytracing%nrayr < 5) then
if (.not. params%raytracing%realtime &
.and. params%raytracing%nrayr < 5) then
params%raytracing%igrad = 0
call log_warning('nrayr < 5 ⇒ optical case only', &
mod="gray_params", proc="set_globals")

2
src/gray_params.sh
View File

@ -11,7 +11,7 @@ sets='antenna equilibrium profiles raytracing ecrh_cd output misc'
antenna='alpha beta power psi chi iox ibeam filenm fghz pos w ri phi'
equilibrium='ssplps ssplf factb sgnb sgni ixp iequil icocos ipsinorm idesc ifreefmt filenm'
profiles='psnbnd sspld factne factte iscal irho iprof filenm'
raytracing='rwmax dst nrayr nrayth nstep igrad idst ipass ipol integrator adaptive_step'
raytracing='rwmax dst nrayr nrayth nstep igrad idst ipass ipol integrator adaptive_step realtime'
ecrh_cd='iwarm ilarm imx ieccd'
output='ipec nrho istpr istpl'
misc='rwall'

47
src/main.f90
View File

@ -32,9 +32,6 @@ program main
if (opts%quiet) opts%verbose = ERROR
call set_log_level(opts%verbose)
! Activate the given output units
call set_active_units(opts%units)
! Load the parameters from file and move to its directory
! (all other filepaths are assumed relative to it)
if (allocated(opts%config_file)) then
@ -65,10 +62,23 @@ program main
! Apply CLI overrides to the parameters
call parse_param_overrides(params)
! Realtime mode
if (params%raytracing%realtime) then
params%ecrh_cd%ieccd = 0 ! Disable current drive compuration
params%ecrh_cd%iwarm = 1 ! Use the weakly relativistic dispersion
params%raytracing%nrayr = 1 ! One ray
params%raytracing%ipass = 1 ! Single pass only
opts%units = [0] ! Disable all output files
call log_message(level=INFO, mod='main', msg='running in realtime mode')
end if
! Copy the parameters into global variables
! exported by the gray_params module
call params_set_globals(params)
! Activate the given output units
call set_active_units(opts%units)
! Read the input data and set the global variables
! of the respective module. Note: order matters.
call init_equilibrium(params, data, err)
@ -93,6 +103,7 @@ program main
end if
if (allocated(opts%sum_filelist)) then
! Combine the output profiles from many individual simulations
call log_message(level=INFO, mod='main', msg='summing profiles')
sum: block
@ -182,23 +193,35 @@ program main
deallocate(opts%params_file)
end block sum
else
! Run the main GRAY routine
call gray_main(params, data, results, err)
end if
print_res: block
print_results: block
character(256) :: msg
write(msg, '(a,g0.5)') 'absorption peak: ρ=', results%rho_peak
call log_message(msg, level=INFO, mod='main')
if (.not. params%raytracing%realtime) then
if (params%ecrh_cd%iwarm > 0) then
write(msg, '(a,g0.3," MW")') 'total absoption: P=', results%pabs
call log_message(msg, level=INFO, mod='main')
write(msg, '(a,g0.3," kA")') 'total current drive: I=', results%icd * 1.0e3_wp_
end if
if (params%ecrh_cd%ieccd > 0) then
write(msg, '(a,g0.3," kA")') 'total current drive: I=', results%icd * 1000
call log_message(msg, level=INFO, mod='main')
end block print_res
end if
end if
end block print_results
! Free memory
call deinit_equilibrium(data%equilibrium)
call deinit_profiles(data%profiles)
call deinit_misc
call deinit_cli_options(opts)
deallocate(results%dpdv, results%jcd)
if (allocated(results%dpdv)) deallocate(results%dpdv, results%jcd)
call close_units
contains
@ -496,7 +519,7 @@ contains
real(wp_), dimension(:, :), pointer :: psjki=>null(), ppabs=>null(), ccci=>null()
real(wp_), dimension(:), pointer :: tau0=>null(), alphaabs0=>null(), &
dids0=>null(), ccci0=>null()
dpds=>null(), dids0=>null(), ccci0=>null()
real(wp_), dimension(:), pointer :: p0jk=>null()
complex(wp_), dimension(:), pointer :: ext=>null(), eyt=>null()
integer, dimension(:), pointer :: iiv=>null()
@ -511,8 +534,8 @@ contains
! Initialise the ray variables (beamtracing)
call init_btr(params%raytracing, yw, ypw, xc, du1, &
gri, ggri, psjki, ppabs, ccci, &
tau0, alphaabs0, dids0, ccci0, &
p0jk, ext, eyt, iiv)
tau0, alphaabs0, dpds, dids0, &
ccci0, p0jk, ext, eyt, iiv)
! Initialise the dispersion module
if (params%ecrh_cd%iwarm > 1) call expinit
@ -554,8 +577,8 @@ contains
! Free memory
call dealloc_surfvec ! for fluxval
call dealloc_beam(yw, ypw, xc, du1, gri, ggri, psjki, ppabs, ccci, &
tau0, alphaabs0, dids0, ccci0, p0jk, ext, eyt, iiv)
call dealloc_beam(yw, ypw, xc, du1, gri, ggri, psjki, ppabs, ccci, tau0, &
alphaabs0, dpds, dids0, ccci0, p0jk, ext, eyt, iiv)
call dealloc_pec
end subroutine sum_profiles

41
src/utils.f90
View File

@ -129,6 +129,47 @@ contains
y=aa*y1+bb*y2
end subroutine intlin
pure function parabola_vertex(x, y) result(v)
! Computes the vertex of a parabola passing by three points
! (x₁, y₁), (x₂, y₂), (x₃, y₃)
!
! The solution is obtained by solving the system
!
! p(x₁) = y₁
! p(x₂) = y₂
! p(x₃) = y₃
!
! where p(x) = ax² + bx + c. The vertex is then
! given by the usual formula (-b/2a, c - b²/4a).
implicit none
! function arguments
real(wp_), dimension(3), intent(in) :: x, y ! x,y of the three points
real(wp_), dimension(2) :: v ! vertex
! local variables
real(wp_) :: denom, a, b, c
denom = (x(1) - x(2))*(x(1) - x(3))*(x(2) - x(3))
a = ( x(3)*(y(2) - y(1)) &
+ x(2)*(y(1) - y(3)) &
+ x(1)*(y(3) - y(2))) / denom
b = ( x(3)**2 * (y(1) - y(2)) &
+ x(2)**2 * (y(3) - y(1)) &
+ x(1)**2 * (y(2) - y(3))) / denom
c = ( x(2)*x(3)*(x(2) - x(3))*y(1) &
+ x(3)*x(1)*(x(3) - x(1))*y(2) &
+ x(1)*x(2)*(x(1) - x(2))*y(3)) / denom
v = [-b/(2*a), c - b**2/(4*a)]
end function
subroutine vmax(x,n,xmax,imx)
implicit none
integer, intent(in) :: n