From 7342566ac051a4bae5a18771ee8e4d7cc3ea2a78 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Tue, 24 May 2022 11:19:30 +0200 Subject: [PATCH] add a switch for realtime mode MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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; --- input/gray.ini | 12 ++++ src/beamdata.f90 | 21 +++--- src/gray_core.f90 | 161 ++++++++++++++++++++++++++++++-------------- src/gray_params.f90 | 6 +- src/gray_params.sh | 2 +- src/main.f90 | 51 ++++++++++---- src/utils.f90 | 41 +++++++++++ 7 files changed, 218 insertions(+), 76 deletions(-) diff --git a/input/gray.ini b/input/gray.ini index 18798d1..84166d5 100644 --- a/input/gray.ini +++ b/input/gray.ini @@ -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 diff --git a/src/beamdata.f90 b/src/beamdata.f90 index a36c84a..0f95230 100644 --- a/src/beamdata.f90 +++ b/src/beamdata.f90 @@ -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) diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 68b3cb9..281da95 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -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,17 +111,35 @@ 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 + if (params%ecrh_cd%iwarm > 1) call expinit - ! Initialise the magsurf_data module - call flux_average ! requires frhotor for dadrhot,dvdrhot + ! 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 + ! Note: 1. needed for computing flux surface averages (dP/dV, etc.) + ! 2. requires frhotor for dadrhot, dvdrhot + call flux_average - ! Initialise the output profiles - call pec_init(params%output%ipec, rhout) - nnd = size(rhop_tab) ! number of radial profile points + ! 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, & @@ -127,34 +147,31 @@ contains 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 ======== - call print_headers(params) + 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 - call print_surfq([1.5_wp_, 2.0_wp_]) + ! print ψ surface for q=1.5 and q=2 on file and log psi,rhot,rhop + call print_surfq([1.5_wp_, 2.0_wp_]) - ! print initial position - write (msg, '("initial position:",3(x,g0.3))') params%antenna%pos - call log_info(msg, mod='gray_core', proc='gray_main') + ! print initial position + write (msg, '("initial position:",3(x,g0.3))') params%antenna%pos + call log_info(msg, mod='gray_core', proc='gray_main') - write (msg, '("initial direction:",2(x,a,"=",g0.2))') & - 'α', params%antenna%alpha, 'β', params%antenna%beta - call log_info(msg, mod='gray_core', proc='gray_main') + write (msg, '("initial direction:",2(x,a,"=",g0.2))') & + 'α', params%antenna%alpha, 'β', params%antenna%beta + call log_info(msg, mod='gray_core', proc='gray_main') - ! print Btot=Bres - ! print ne, Te, q, Jphi versus psi, rhop, rhot - call print_bres(bres) - call print_prof(params%profiles) - call print_maps(bres, xgcn, & - norm2(params%antenna%pos(1:2)) * 0.01_wp_, & - sin(params%antenna%beta*degree)) + ! print Btot=Bres + ! print ne, Te, q, Jphi versus psi, rhop, rhot + call print_bres(bres) + call print_prof(params%profiles) + 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 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 diff --git a/src/gray_params.f90 b/src/gray_params.f90 index c29f1f5..41c057e 100644 --- a/src/gray_params.f90 +++ b/src/gray_params.f90 @@ -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") diff --git a/src/gray_params.sh b/src/gray_params.sh index 76da5dc..36c0b77 100644 --- a/src/gray_params.sh +++ b/src/gray_params.sh @@ -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' diff --git a/src/main.f90 b/src/main.f90 index c4f9730..17a64bd 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -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.3," MW")') 'total absoption: P=', results%pabs + + write(msg, '(a,g0.5)') 'absorption peak: ρ=', results%rho_peak call log_message(msg, level=INFO, mod='main') - write(msg, '(a,g0.3," kA")') 'total current drive: I=', results%icd * 1.0e3_wp_ - call log_message(msg, level=INFO, mod='main') - end block print_res + + 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') + 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 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 diff --git a/src/utils.f90 b/src/utils.f90 index a2cc142..0863b5e 100644 --- a/src/utils.f90 +++ b/src/utils.f90 @@ -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