gray-jintrac gray-externals.f: fixed fghz error in read_data, updated read_beams
This commit is contained in:
Daniele Micheletti
parent
b9b687be88
commit
72b3c5f511
@ -679,21 +679,14 @@ c
|
|||||||
end if
|
end if
|
||||||
c
|
c
|
||||||
fhz=fghz*1.0d9
|
fhz=fghz*1.0d9
|
||||||
ak0=2.0d9*pi*fghz/vc
|
|
||||||
akinv=1.0d0/ak0
|
|
||||||
c
|
|
||||||
bresg=2.0d0*pi*fhz*me*vc/qe
|
|
||||||
bres=bresg*1.0d-4
|
|
||||||
c
|
|
||||||
c xg=xgcn*dens19
|
|
||||||
c
|
|
||||||
xgcn=1.0d-5*qe**2/(pi*me*fghz**2)
|
|
||||||
c
|
c
|
||||||
c read data for beam from file if ibeam>0
|
c read data for beam from file if ibeam>0
|
||||||
c
|
c
|
||||||
if(ibeam.gt.0) then
|
if(ibeam.gt.0) then
|
||||||
call read_beams(beamid,iox)
|
call read_beams(beamid,iox)
|
||||||
|
ak0=2.0d9*pi*fghz/vc
|
||||||
else
|
else
|
||||||
|
ak0=2.0d9*pi*fghz/vc
|
||||||
zrcsi=0.5d0*ak0*w0csi**2
|
zrcsi=0.5d0*ak0*w0csi**2
|
||||||
zreta=0.5d0*ak0*w0eta**2
|
zreta=0.5d0*ak0*w0eta**2
|
||||||
if(igrad.gt.0) then
|
if(igrad.gt.0) then
|
||||||
@ -716,6 +709,15 @@ c
|
|||||||
end if
|
end if
|
||||||
end if
|
end if
|
||||||
c
|
c
|
||||||
|
akinv=1.0d0/ak0
|
||||||
|
c
|
||||||
|
bresg=2.0d0*pi*fhz*me*vc/qe
|
||||||
|
bres=bresg*1.0d-4
|
||||||
|
c
|
||||||
|
c xg=xgcn*dens19
|
||||||
|
c
|
||||||
|
xgcn=1.0d-5*qe**2/(pi*me*fghz**2)
|
||||||
|
c
|
||||||
sox=-1.0d0
|
sox=-1.0d0
|
||||||
if(iox.eq.2) sox=1.0d0
|
if(iox.eq.2) sox=1.0d0
|
||||||
c
|
c
|
||||||
@ -894,34 +896,31 @@ c
|
|||||||
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
|
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
|
||||||
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
|
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
|
||||||
subroutine read_beams(beamid,iox)
|
subroutine read_beams(beamid,iox)
|
||||||
|
use const_and_precisions, only : pi
|
||||||
use reflections, only : inside
|
use reflections, only : inside
|
||||||
implicit real*8(a-h,o-z)
|
implicit real*8(a-h,o-z)
|
||||||
c
|
c
|
||||||
character*255 filenmeqq,filenmprf,filenmbm
|
character*255 filenmeqq,filenmprf,filenmbm
|
||||||
character*20 beamname,beamnameskip
|
character*20 beamname
|
||||||
c
|
c
|
||||||
integer beamid, nisteer, fdeg, jumprow, nbeam, nalpha, nbeta,
|
integer beamid, nisteer, fdeg, jumprow, nbeam, nalpha, nbeta,
|
||||||
. iopt, incheck, nxcoord, nycoord, nxest, nyest, lwrk, kwrk,
|
. iopt, incheck, nxcoord, nycoord, nxest, nyest, lwrk, kwrk,
|
||||||
. nxwaist1, nywaist1, nxwaist2, nywaist2, nxrci1, nyrci1, nxrci2,
|
. nx, ny, kx, ky, ii, nfbeam, iox, npolyg, nmax, lwrk2,
|
||||||
. nyrci2, nxphi1, nyphi1, nxphi2, nyphi2, nxx0, nyx0, nxy0, nyy0,
|
. in, lwrk_ev, kwrk_ev
|
||||||
. nxz0, nyz0, kx, ky, ii, nfbeam, iox, npolyg, nmax, lwrk2, in
|
integer, DIMENSION(:), ALLOCATABLE :: iwrk, iwrk_ev
|
||||||
integer, DIMENSION(:), ALLOCATABLE :: iwrk
|
|
||||||
c
|
c
|
||||||
real*8 alphast, betast, x00, y00, z00, waist1, waist2, dvir1,
|
real*8 alphast, betast, x00, y00, z00, waist1, waist2, dvir1,
|
||||||
. dvir2, delta, phi1, phi, zr1, zr2, ako, w1, w2, rci1, rci2, fghz,
|
. dvir2, delta, phi1, phi2, zr1, zr2, w1, w2, rci1, rci2, fghz,
|
||||||
. fhz, wcso, weta, rcicsi, rcieta, phiw, phir, xcoord0, ycoord0,
|
. fhz, wcsi, weta, rcicsi, rcieta, phiw, phir, xcoord0, ycoord0,
|
||||||
. alpha0, beta0, fp, minx, maxx, miny, maxy, eps, dmin
|
. alpha0, beta0, fp, minx, maxx, miny, maxy, eps, dmin, vc
|
||||||
real*8, DIMENSION(:), ALLOCATABLE :: x00v, y00v, z00v, alphastv,
|
real*8, DIMENSION(:), ALLOCATABLE :: x00v, y00v, z00v, alphastv,
|
||||||
. betastv, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v, xcoord,
|
. betastv, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v, xcoord,
|
||||||
. ycoord, ycoord2, wrk, txwaist1, tywaist1, txwaist2, tywaist2,
|
. ycoord, wrk, wrk_ev, tx, ty, c, w, wrk2, xpolyg, ypolyg, xpolygA,
|
||||||
. txrci1, tyrci1, txrci2, tyrci2, txphi1, typhi1, txphi2, typhi2,
|
. ypolygA, xpolygB, ypolygB, xpolygC, ypolygC, xpolygD, ypolygD
|
||||||
. txx0, tyx0, txy0, tyy0, txz0, tyz0, txycoord, cycoord, cwaist1,
|
|
||||||
. cwaist2, crci1, crci2, cphi1,cphi2, cx0, cy0, cz0, w, wrk2,
|
|
||||||
. xpolyg, ypolyg, xpolygA, ypolygA, xpolygB, ypolygB, xpolygC,
|
|
||||||
. ypolygC, xpolygD, ypolygD
|
|
||||||
real*8 xvert(4), yvert(4)
|
real*8 xvert(4), yvert(4)
|
||||||
c
|
c
|
||||||
parameter(kspl=1,sspl=0.01)
|
parameter(kspl=1,sspl=0.01d0)
|
||||||
|
parameter(vc=2.9979d+10)
|
||||||
c
|
c
|
||||||
common/ibeam/ibeam
|
common/ibeam/ibeam
|
||||||
common/filesn/filenmeqq,filenmprf,filenmbm
|
common/filesn/filenmeqq,filenmprf,filenmbm
|
||||||
@ -941,6 +940,7 @@ c
|
|||||||
c
|
c
|
||||||
c#######################################################################################
|
c#######################################################################################
|
||||||
if(ibeam.eq.1) then
|
if(ibeam.eq.1) then
|
||||||
|
ak0=2.0d9*pi*fghz/vc
|
||||||
read(nfbeam,*) nisteer
|
read(nfbeam,*) nisteer
|
||||||
naplha = nisteer
|
naplha = nisteer
|
||||||
nbeta = 1
|
nbeta = 1
|
||||||
@ -985,11 +985,11 @@ c===============================================================================
|
|||||||
c # of beams
|
c # of beams
|
||||||
read(nfbeam,*) nbeam
|
read(nfbeam,*) nbeam
|
||||||
c
|
c
|
||||||
c unused beams' data
|
c unused beam data
|
||||||
jumprow=0
|
jumprow=0
|
||||||
c c====================================================================================
|
c c====================================================================================
|
||||||
do i=1,beamid-1
|
do i=1,beamid-1
|
||||||
read(nfbeam,*) beamnameskip, iox, fghz, nalpha, nbeta
|
read(nfbeam,*) beamname, iox, fghz, nalpha, nbeta
|
||||||
jumprow = jumprow+nalpha*nbeta
|
jumprow = jumprow+nalpha*nbeta
|
||||||
end do
|
end do
|
||||||
c c====================================================================================
|
c c====================================================================================
|
||||||
@ -999,13 +999,9 @@ c beam of interest
|
|||||||
fhz=fghz*1.0d9
|
fhz=fghz*1.0d9
|
||||||
c
|
c
|
||||||
c c====================================================================================
|
c c====================================================================================
|
||||||
c unused beams' data grids
|
c unused beam data grids
|
||||||
do i=1,(nbeam - beamid)
|
do i=1,(jumprow + (nbeam - beamid))
|
||||||
read(nfbeam,*) beamnameskip
|
read(nfbeam,*)
|
||||||
end do
|
|
||||||
do i=1,jumprow
|
|
||||||
read(nfbeam,*) alphast,betast,x00,y00,z00,
|
|
||||||
. waist1,waist2,rci1,rci2,phi1,phi2
|
|
||||||
end do
|
end do
|
||||||
c c====================================================================================
|
c c====================================================================================
|
||||||
c
|
c
|
||||||
@ -1095,63 +1091,10 @@ c
|
|||||||
iopt = 0
|
iopt = 0
|
||||||
incheck = 0
|
incheck = 0
|
||||||
c
|
c
|
||||||
|
c#######################################################################################
|
||||||
|
c POINT POSITION CHECK
|
||||||
c#######################################################################################
|
c#######################################################################################
|
||||||
if(fdeg.ne.0) then
|
if(fdeg.ne.0) then
|
||||||
nxest = kx + nxcoord + 1
|
|
||||||
lwrk = (nxcoord*(kx+1)+nxest*(7+3*kx))
|
|
||||||
kwrk = nxest
|
|
||||||
allocate(cycoord(nxest), txycoord(nxest), cwaist1(nxest),
|
|
||||||
. txwaist1(nxest), cwaist2(nxest), txwaist2(nxest),
|
|
||||||
. crci1(nxest), txrci1(nxest), crci2(nxest), txrci2(nxest),
|
|
||||||
. cphi1(nxest), txphi1(nxest), cphi2(nxest), txphi2(nxest),
|
|
||||||
. cx0(nxest), txx0(nxest), cy0(nxest), txy0(nxest),
|
|
||||||
. cz0(nxest), txz0(nxest), w(nxcoord), wrk(lwrk), iwrk(kwrk))
|
|
||||||
c
|
|
||||||
do i=1,nxcoord
|
|
||||||
w(i) = 1.d0
|
|
||||||
end do
|
|
||||||
c
|
|
||||||
c 2D interpolation
|
|
||||||
call dierckx_curfit(iopt,nxcoord,xcoord,ycoord,w,
|
|
||||||
. xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxycoord,
|
|
||||||
. txycoord,cycoord,fp,wrk,lwrk,iwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_curfit(iopt,nxcoord,xcoord,waist1v,w,
|
|
||||||
. xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxwaist1,
|
|
||||||
. txwaist1,cwaist1,fp,wrk,lwrk,iwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_curfit(iopt,nxcoord,xcoord,waist2v,w,
|
|
||||||
. xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxwaist2,
|
|
||||||
. txwaist2,cwaist2,fp,wrk,lwrk,iwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_curfit(iopt,nxcoord,xcoord,rci1v,w,
|
|
||||||
. xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxrci1,
|
|
||||||
. txrci1,crci1,fp,wrk,lwrk,iwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_curfit(iopt,nxcoord,xcoord,rci2v,w,
|
|
||||||
. xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxrci2,
|
|
||||||
. txrci2,crci2,fp,wrk,lwrk,iwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_curfit(iopt,nxcoord,xcoord,phi1v,w,
|
|
||||||
. xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxphi1,
|
|
||||||
. txphi1,cphi1,fp,wrk,lwrk,iwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_curfit(iopt,nxcoord,xcoord,phi2v,w,
|
|
||||||
. xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxphi2,
|
|
||||||
. txphi2,cphi2,fp,wrk,lwrk,iwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_curfit(iopt,nxcoord,xcoord,x00v,w,
|
|
||||||
. xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxx0,
|
|
||||||
. txx0,cx0,fp,wrk,lwrk,iwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_curfit(iopt,nxcoord,xcoord,y00v,w,
|
|
||||||
. xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxy0,
|
|
||||||
. txy0,cy0,fp,wrk,lwrk,iwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_curfit(iopt,nxcoord,xcoord,z00v,w,
|
|
||||||
. xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxz0,
|
|
||||||
. txz0,cz0,fp,wrk,lwrk,iwrk,ier)
|
|
||||||
c
|
|
||||||
c check if xcoord0 is out of table range
|
c check if xcoord0 is out of table range
|
||||||
c incheck = 1 inside / 0 outside
|
c incheck = 1 inside / 0 outside
|
||||||
if(xcoord0.gt.xcoord(1).and.xcoord0.lt.xcoord(nisteer))
|
if(xcoord0.gt.xcoord(1).and.xcoord0.lt.xcoord(nisteer))
|
||||||
@ -1160,83 +1103,11 @@ c c============================================================================
|
|||||||
else
|
else
|
||||||
c c====================================================================================
|
c c====================================================================================
|
||||||
npolyg = 2*(nxcoord+nycoord-2)
|
npolyg = 2*(nxcoord+nycoord-2)
|
||||||
minx = minval(xcoord)
|
allocate(xpolyg(npolyg), ypolyg(npolyg))
|
||||||
maxx = maxval(xcoord)
|
|
||||||
miny = minval(ycoord)
|
|
||||||
maxy = maxval(ycoord)
|
|
||||||
nxest = kx + 1 + sqrt(nisteer/2.)
|
|
||||||
nyest = ky + 1 + sqrt(nisteer/2.)
|
|
||||||
nmax = max(nxest,nyest)
|
|
||||||
eps = 10.**(-8)
|
|
||||||
lwrk = (nmax-2)**2*(7*nmax-2)+18*nmax+8*nisteer-19
|
|
||||||
lwrk2 = (nmax-2)**2*(4*nmax-1)+4*nmax-2
|
|
||||||
kwrk = nisteer+(nmax-3)*(nmax-3)
|
|
||||||
allocate(cwaist1(nxest*nyest), txwaist1(nmax), tywaist1(nmax),
|
|
||||||
. cwaist2(nxest*nyest), txwaist2(nmax), tywaist2(nmax),
|
|
||||||
. crci1(nxest*nyest), txrci1(nmax), tyrci1(nmax),
|
|
||||||
. crci2(nxest*nyest), txrci2(nmax), tyrci2(nmax),
|
|
||||||
. cphi1(nxest*nyest), txphi1(nmax), typhi1(nmax),
|
|
||||||
. cphi2(nxest*nyest), txphi2(nmax), typhi2(nmax),
|
|
||||||
. cx0(nxest*nyest), txx0(nmax), tyx0(nmax),
|
|
||||||
. cy0(nxest*nyest), txy0(nmax), tyy0(nmax),
|
|
||||||
. cz0(nxest*nyest), txz0(nmax), tyz0(nmax),
|
|
||||||
. wrk(lwrk), wrk2(lwrk2), iwrk(kwrk),
|
|
||||||
. xpolyg(npolyg), ypolyg(npolyg), w(nisteer))
|
|
||||||
c
|
|
||||||
do i=1,nisteer
|
|
||||||
w(i) = 1.d0
|
|
||||||
end do
|
|
||||||
c
|
|
||||||
c 3D interpolation
|
|
||||||
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,waist1v,w,
|
|
||||||
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
|
||||||
. nxwaist1,txwaist1,nywaist1,tywaist1,cwaist1,
|
|
||||||
. fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,waist2v,w,
|
|
||||||
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
|
||||||
. nxwaist2,txwaist2,nywaist2,tywaist2,cwaist2,
|
|
||||||
. fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,rci1v,w,
|
|
||||||
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
|
||||||
. nxrci1,txrci1,nyrci1,tyrci1,crci1,
|
|
||||||
. fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,rci2v,w,
|
|
||||||
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
|
||||||
. nxrci2,txrci2,nyrci2,tyrci2,crci2,
|
|
||||||
. fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,phi1v,w,
|
|
||||||
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
|
||||||
. nxphi1,txphi1,nyphi1,typhi1,cphi1,
|
|
||||||
. fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,phi2v,w,
|
|
||||||
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
|
||||||
. nxphi2,txphi2,nyphi2,typhi2,cphi2,
|
|
||||||
. fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,x00v,w,
|
|
||||||
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
|
||||||
. nxx0,txx0,nyx0,tyx0,cx0,
|
|
||||||
. fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,y00v,w,
|
|
||||||
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
|
||||||
. nxy0,txy0,nyy0,tyy0,cy0,
|
|
||||||
. fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,z00v,w,
|
|
||||||
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
|
||||||
. nxz0,txz0,nyz0,tyz0,cz0,
|
|
||||||
. fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
|
||||||
c data range polygon
|
c data range polygon
|
||||||
xpolyg(1:nxcoord) = xcoord(1:nxcoord)
|
xpolyg(1:nxcoord) = xcoord(1:nxcoord)
|
||||||
ypolyg(1:nxcoord) = ycoord(1:nxcoord)
|
ypolyg(1:nxcoord) = ycoord(1:nxcoord)
|
||||||
c
|
c
|
||||||
c c====================================================================================
|
|
||||||
do i=1,nycoord-2
|
do i=1,nycoord-2
|
||||||
xpolyg(nxcoord+i) = xcoord((i+1)*nxcoord)
|
xpolyg(nxcoord+i) = xcoord((i+1)*nxcoord)
|
||||||
xpolyg(2*nxcoord+nycoord-2+i) =
|
xpolyg(2*nxcoord+nycoord-2+i) =
|
||||||
@ -1257,41 +1128,14 @@ c incheck = 1 inside / 0 outside
|
|||||||
if(inside(xpolyg,ypolyg,npolyg,xcoord0,ycoord0))
|
if(inside(xpolyg,ypolyg,npolyg,xcoord0,ycoord0))
|
||||||
. incheck = 1
|
. incheck = 1
|
||||||
end if
|
end if
|
||||||
deallocate(wrk,iwrk)
|
|
||||||
c#######################################################################################
|
c#######################################################################################
|
||||||
c
|
c
|
||||||
c
|
c#######################################################################################
|
||||||
|
c INTERPOLATION
|
||||||
c#######################################################################################
|
c#######################################################################################
|
||||||
if(fdeg.ne.0) then
|
if(fdeg.ne.0) then
|
||||||
c c====================================================================================
|
c c====================================================================================
|
||||||
if(incheck.eq.1) then
|
if(incheck.eq.0) then
|
||||||
call dierckx_splev(txycoord,nxycoord,cycoord,kx,xcoord0,
|
|
||||||
. ycoord0,1,ier)
|
|
||||||
c
|
|
||||||
call dierckx_splev(txwaist1,nxwaist1,cwaist1,kx,xcoord0,wcsi,
|
|
||||||
. 1,ier)
|
|
||||||
c
|
|
||||||
call dierckx_splev(txwaist2,nxwaist2,cwaist2,kx,xcoord0,weta,
|
|
||||||
. 1,ier)
|
|
||||||
c
|
|
||||||
call dierckx_splev(txrci1,nxrci1,crci1,kx,xcoord0,rcicsi,
|
|
||||||
. 1,ier)
|
|
||||||
c
|
|
||||||
call dierckx_splev(txrci2,nxrci2,crci2,kx,xcoord0,rcieta,
|
|
||||||
. 1,ier)
|
|
||||||
c
|
|
||||||
call dierckx_splev(txphi1,nxphi1,cphi1,kx,xcoord0,phiw,1,ier)
|
|
||||||
c
|
|
||||||
call dierckx_splev(txphi2,nxphi2,cphi2,kx,xcoord0,phir,1,ier)
|
|
||||||
c
|
|
||||||
call dierckx_splev(txx0,nxx0,cx0,kx,xcoord0,x00,1,ier)
|
|
||||||
c
|
|
||||||
call dierckx_splev(txy0,nxy0,cy0,kx,xcoord0,y00,1,ier)
|
|
||||||
c
|
|
||||||
call dierckx_splev(txz0,nxz0,cz0,kx,xcoord0,z00,1,ier)
|
|
||||||
c c----------------------------------------------------------------------------------
|
|
||||||
else
|
|
||||||
c c----------------------------------------------------------------------------------
|
|
||||||
write(*,*) ' alpha0 or beta0 outside table range !!!'
|
write(*,*) ' alpha0 or beta0 outside table range !!!'
|
||||||
if(xcoord0.ge.xcoord(nisteer)) ii=nisteer
|
if(xcoord0.ge.xcoord(nisteer)) ii=nisteer
|
||||||
if(xcoord0.le.xcoord(1)) ii=1
|
if(xcoord0.le.xcoord(1)) ii=1
|
||||||
@ -1307,10 +1151,73 @@ c
|
|||||||
rcieta=rci2v(ii)
|
rcieta=rci2v(ii)
|
||||||
phiw=phi1v(ii)
|
phiw=phi1v(ii)
|
||||||
phir=phi2v(ii)
|
phir=phi2v(ii)
|
||||||
|
c c----------------------------------------------------------------------------------
|
||||||
|
else
|
||||||
|
c c----------------------------------------------------------------------------------
|
||||||
|
nxest = kx + nxcoord + 1
|
||||||
|
lwrk = (nxcoord*(kx+1)+nxest*(7+3*kx))
|
||||||
|
kwrk = nxest
|
||||||
|
allocate(c(nxest),tx(nxest),w(nxcoord),wrk(lwrk),iwrk(kwrk))
|
||||||
|
c
|
||||||
|
do i=1,nxcoord
|
||||||
|
w(i) = 1.d0
|
||||||
|
end do
|
||||||
|
c
|
||||||
|
c coefficients calculation & evaluation in xcoord0
|
||||||
|
call dierckx_curfit(iopt,nxcoord,xcoord,ycoord,w,xcoord(1),
|
||||||
|
. xcoord(nxcoord),kx,sspl,nxest,nx,tx,c,fp,wrk,lwrk,iwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_splev(tx,nx,c,kx,xcoord0,ycoord0,1,ier)
|
||||||
|
c
|
||||||
|
call dierckx_curfit(iopt,nxcoord,xcoord,waist1v,w,xcoord(1),
|
||||||
|
. xcoord(nxcoord),kx,sspl,nxest,nx,tx,c,fp,wrk,lwrk,iwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_splev(tx,nx,c,kx,xcoord0,wcsi,1,ier)
|
||||||
|
c
|
||||||
|
call dierckx_curfit(iopt,nxcoord,xcoord,waist2v,w,xcoord(1),
|
||||||
|
. xcoord(nxcoord),kx,sspl,nxest,nx,tx,c,fp,wrk,lwrk,iwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_splev(tx,nx,c,kx,xcoord0,weta,1,ier)
|
||||||
|
c
|
||||||
|
call dierckx_curfit(iopt,nxcoord,xcoord,rci1v,w,xcoord(1),
|
||||||
|
. xcoord(nxcoord),kx,sspl,nxest,nx,tx,c,fp,wrk,lwrk,iwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_splev(tx,nx,c,kx,xcoord0,rcicsi,1,ier)
|
||||||
|
c
|
||||||
|
call dierckx_curfit(iopt,nxcoord,xcoord,rci2v,w,xcoord(1),
|
||||||
|
. xcoord(nxcoord),kx,sspl,nxest,nx,tx,c,fp,wrk,lwrk,iwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_splev(tx,nx,c,kx,xcoord0,rcieta,1,ier)
|
||||||
|
c
|
||||||
|
call dierckx_curfit(iopt,nxcoord,xcoord,phi1v,w,xcoord(1),
|
||||||
|
. xcoord(nxcoord),kx,sspl,nxest,nx,tx,c,fp,wrk,lwrk,iwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_splev(tx,nx,c,kx,xcoord0,phiw,1,ier)
|
||||||
|
c
|
||||||
|
call dierckx_curfit(iopt,nxcoord,xcoord,phi2v,w,xcoord(1),
|
||||||
|
. xcoord(nxcoord),kx,sspl,nxest,nx,tx,c,fp,wrk,lwrk,iwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_splev(tx,nx,c,kx,xcoord0,phir,1,ier)
|
||||||
|
c
|
||||||
|
call dierckx_curfit(iopt,nxcoord,xcoord,x00v,w,xcoord(1),
|
||||||
|
. xcoord(nxcoord),kx,sspl,nxest,nx,tx,c,fp,wrk,lwrk,iwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_splev(tx,nx,c,kx,xcoord0,x00,1,ier)
|
||||||
|
c
|
||||||
|
call dierckx_curfit(iopt,nxcoord,xcoord,y00v,w,xcoord(1),
|
||||||
|
. xcoord(nxcoord),kx,sspl,nxest,nx,tx,c,fp,wrk,lwrk,iwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_splev(tx,nx,c,kx,xcoord0,y00,1,ier)
|
||||||
|
c
|
||||||
|
call dierckx_curfit(iopt,nxcoord,xcoord,z00v,w,xcoord(1),
|
||||||
|
. xcoord(nxcoord),kx,sspl,nxest,nx,tx,c,fp,wrk,lwrk,iwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_splev(tx,nx,c,kx,xcoord0,z00,1,ier)
|
||||||
end if
|
end if
|
||||||
c c====================================================================================
|
c c====================================================================================
|
||||||
else
|
else
|
||||||
c c====================================================================================
|
c c====================================================================================
|
||||||
|
c the nearest grid point for a (xcoord0,ycoord0) external point is found
|
||||||
if(incheck.eq.0) then
|
if(incheck.eq.0) then
|
||||||
allocate(xpolygA(nxcoord), ypolygA(nxcoord), xpolygC(nxcoord),
|
allocate(xpolygA(nxcoord), ypolygA(nxcoord), xpolygC(nxcoord),
|
||||||
. ypolygC(nxcoord), xpolygB(nycoord), ypolygB(nycoord),
|
. ypolygC(nxcoord), xpolygB(nycoord), ypolygB(nycoord),
|
||||||
@ -1376,15 +1283,15 @@ c
|
|||||||
c c----------------------------------------------------------------------------------
|
c c----------------------------------------------------------------------------------
|
||||||
c
|
c
|
||||||
c c----------------------------------------------------------------------------------
|
c c----------------------------------------------------------------------------------
|
||||||
c (xcoord0,ycoord0) is set to its nearest point on (alpha, beta) grid border
|
c (xcoord0,ycoord0) is set to a point on (alpha, beta) grid border
|
||||||
c depending on the region
|
c depending on the region
|
||||||
c 1: xcoord0 unchanged, ycoord0 moved on side A
|
c 1: xcoord0 unchanged, ycoord0 moved on side A
|
||||||
c 3: xcoord0 moved on side B, ycoord0 unchanged
|
c 3: xcoord0 moved on side B, ycoord0 unchanged
|
||||||
c 5: xcoord0 unchanged, ycoord0 moved on side C
|
c 5: xcoord0 unchanged, ycoord0 moved on side C
|
||||||
c 7: xcoord0 moved on side D, ycoord0 unchanged
|
c 7: xcoord0 moved on side D, ycoord0 unchanged
|
||||||
c 2,4,6,8: (xcoord0,ycoord0) set to nearest vertex coordinates
|
c 2,4,6,8: (xcoord0,ycoord0) set to nearest vertex coordinates
|
||||||
c in 1,3,5,7 incheck is set back to 1 to evaluate x00,y00,z00,waist,rci,phi in
|
c in 1,3,5,7 incheck is set back to 1 to evaluate x00,y00,z00,waist,rci,phi
|
||||||
c new (xcoord0,ycoord0)
|
c at the new (xcoord0,ycoord0) coordinates
|
||||||
c in 2,4,6,8 incheck remains 0 and x00,y00,z00,waist,rci,phi values at the
|
c in 2,4,6,8 incheck remains 0 and x00,y00,z00,waist,rci,phi values at the
|
||||||
c (xcoord0,ycoord0) vertex are used
|
c (xcoord0,ycoord0) vertex are used
|
||||||
alpha0 = xcoord0
|
alpha0 = xcoord0
|
||||||
@ -1438,6 +1345,9 @@ c xcoord0, ycoord0 set
|
|||||||
ycoord0 = yvert(1)
|
ycoord0 = yvert(1)
|
||||||
ii = 1
|
ii = 1
|
||||||
END SELECT
|
END SELECT
|
||||||
|
c
|
||||||
|
deallocate(xpolygA, ypolygA, xpolygC,ypolygC, xpolygB,
|
||||||
|
. ypolygB, xpolygD, ypolygD)
|
||||||
c c----------------------------------------------------------------------------------
|
c c----------------------------------------------------------------------------------
|
||||||
c
|
c
|
||||||
c print new (alpha0, beta0)
|
c print new (alpha0, beta0)
|
||||||
@ -1447,49 +1357,11 @@ c print new (alpha0, beta0)
|
|||||||
write(*,10) ' new (alpha0,beta0) = (', xcoord0,
|
write(*,10) ' new (alpha0,beta0) = (', xcoord0,
|
||||||
. ',', ycoord0, ')'
|
. ',', ycoord0, ')'
|
||||||
10 format(a27,f10.5,a1,f10.5,a1)
|
10 format(a27,f10.5,a1,f10.5,a1)
|
||||||
deallocate(xpolygA, ypolygA, xpolygC, ypolygC,
|
|
||||||
. xpolygB, ypolygB, xpolygD, ypolygD)
|
|
||||||
end if
|
end if
|
||||||
c c====================================================================================
|
c c====================================================================================
|
||||||
c
|
c
|
||||||
c c====================================================================================
|
c c====================================================================================
|
||||||
if(incheck.eq.1) then
|
if(incheck.eq.0) then
|
||||||
lwrk = 2*(kx+ky+2)
|
|
||||||
kwrk = 4
|
|
||||||
allocate(wrk(lwrk),iwrk(kwrk))
|
|
||||||
c
|
|
||||||
call dierckx_bispev(txwaist1,nxwaist1,tywaist1,nywaist1,
|
|
||||||
. cwaist1,kx,ky,(/xcoord0/),1,(/ycoord0/),1,wcsi,
|
|
||||||
. wrk,lwrk,iwrk,kwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_bispev(txwaist2,nxwaist2,tywaist2,nywaist2,
|
|
||||||
. cwaist2,kx,ky,(/xcoord0/),1,(/ycoord0/),1,weta,
|
|
||||||
. wrk,lwrk,iwrk,kwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_bispev(txrci1,nxrci1,tyrci1,nyrci1,crci1,kx,
|
|
||||||
. ky,(/xcoord0/),1,(/ycoord0/),1,rcicsi,wrk,lwrk,iwrk,kwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_bispev(txrci2,nxrci2,tyrci2,nyrci2,crci2,kx,
|
|
||||||
. ky,(/xcoord0/),1,(/ycoord0/),1,rcieta,wrk,lwrk,iwrk,kwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_bispev(txphi1,nxphi1,typhi1,nyphi1,cphi1,kx,
|
|
||||||
. ky,(/xcoord0/),1,(/ycoord0/),1,phiw,wrk,lwrk,iwrk,kwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_bispev(txphi2,nxphi2,typhi2,nyphi2,cphi2,kx,
|
|
||||||
. ky,(/xcoord0/),1,(/ycoord0/),1,phir,wrk,lwrk,iwrk,kwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_bispev(txx0,nxx0,tyx0,nyx0,cx0,kx,ky,
|
|
||||||
. (/xcoord0/),1,(/ycoord0/),1,x00,wrk,lwrk,iwrk,kwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_bispev(txy0,nxy0,tyy0,nyy0,cy0,kx,ky,
|
|
||||||
. (/xcoord0/),1,(/ycoord0/),1,y00,wrk,lwrk,iwrk,kwrk,ier)
|
|
||||||
c
|
|
||||||
call dierckx_bispev(txz0,nxz0,tyz0,nyz0,cz0,kx,ky,
|
|
||||||
. (/xcoord0/),1,(/ycoord0/),1,z00,wrk,lwrk,iwrk,kwrk,ier)
|
|
||||||
deallocate(wrk,iwrk)
|
|
||||||
c c----------------------------------------------------------------------------------
|
|
||||||
else
|
|
||||||
c c----------------------------------------------------------------------------------
|
|
||||||
x00=x00v(ii)
|
x00=x00v(ii)
|
||||||
y00=y00v(ii)
|
y00=y00v(ii)
|
||||||
z00=z00v(ii)
|
z00=z00v(ii)
|
||||||
@ -1499,22 +1371,109 @@ c c--------------------------------------------------------------------------
|
|||||||
rcieta=rci2v(ii)
|
rcieta=rci2v(ii)
|
||||||
phiw=phi1v(ii)
|
phiw=phi1v(ii)
|
||||||
phir=phi2v(ii)
|
phir=phi2v(ii)
|
||||||
|
c c----------------------------------------------------------------------------------
|
||||||
|
else
|
||||||
|
c c----------------------------------------------------------------------------------
|
||||||
|
minx = minval(xcoord)
|
||||||
|
maxx = maxval(xcoord)
|
||||||
|
miny = minval(ycoord)
|
||||||
|
maxy = maxval(ycoord)
|
||||||
|
nxest = kx + 1 + sqrt(nisteer/2.)
|
||||||
|
nyest = ky + 1 + sqrt(nisteer/2.)
|
||||||
|
nmax = max(nxest,nyest)
|
||||||
|
eps = 10.**(-8)
|
||||||
|
lwrk = (nmax-2)**2*(7*nmax-2)+18*nmax+8*nisteer-19
|
||||||
|
lwrk2 = (nmax-2)**2*(4*nmax-1)+4*nmax-2
|
||||||
|
kwrk = nisteer+(nmax-3)*(nmax-3)
|
||||||
|
lwrk_ev = 2*(kx+ky+2)
|
||||||
|
kwrk_ev = 4
|
||||||
|
c
|
||||||
|
allocate(c(nxest*nyest), tx(nmax), ty(nmax),
|
||||||
|
. wrk(lwrk), wrk2(lwrk2), iwrk(kwrk), w(nisteer),
|
||||||
|
. wrk_ev(lwrk_ev),iwrk_ev(kwrk_ev))
|
||||||
|
c
|
||||||
|
do i=1,nisteer
|
||||||
|
w(i) = 1.d0
|
||||||
|
end do
|
||||||
|
c
|
||||||
|
c coefficients calculation & evaluation in (xcoord0,ycoord0)
|
||||||
|
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,waist1v,w,
|
||||||
|
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
||||||
|
. nx,tx,ny,ty,c,fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_bispev(tx,nx,ty,ny,c,kx,ky,(/xcoord0/),1,
|
||||||
|
. (/ycoord0/),1,wcsi,wrk_ev,lwrk_ev,iwrk_ev,kwrk_ev,ier)
|
||||||
|
c
|
||||||
|
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,waist2v,w,
|
||||||
|
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
||||||
|
. nx,tx,ny,ty,c,fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_bispev(tx,nx,ty,ny,c,kx,ky,(/xcoord0/),1,
|
||||||
|
. (/ycoord0/),1,weta,wrk_ev,lwrk_ev,iwrk_ev,kwrk_ev,ier)
|
||||||
|
c
|
||||||
|
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,rci1v,w,
|
||||||
|
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
||||||
|
. nx,tx,ny,ty,c,fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_bispev(tx,nx,ty,ny,c,kx,ky,(/xcoord0/),1,
|
||||||
|
. (/ycoord0/),1,rcicsi,wrk_ev,lwrk_ev,iwrk_ev,kwrk_ev,ier)
|
||||||
|
c
|
||||||
|
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,rci2v,w,
|
||||||
|
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
||||||
|
. nx,tx,ny,ty,c,fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_bispev(tx,nx,ty,ny,c,kx,ky,(/xcoord0/),1,
|
||||||
|
. (/ycoord0/),1,rcieta,wrk_ev,lwrk_ev,iwrk_ev,kwrk_ev,ier)
|
||||||
|
c
|
||||||
|
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,phi1v,w,
|
||||||
|
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
||||||
|
. nx,tx,ny,ty,c,fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_bispev(tx,nx,ty,ny,c,kx,ky,(/xcoord0/),1,
|
||||||
|
. (/ycoord0/),1,phiw,wrk_ev,lwrk_ev,iwrk_ev,kwrk_ev,ier)
|
||||||
|
c
|
||||||
|
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,phi2v,w,
|
||||||
|
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
||||||
|
. nx,tx,ny,ty,c,fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_bispev(tx,nx,ty,ny,c,kx,ky,(/xcoord0/),1,
|
||||||
|
. (/ycoord0/),1,phir,wrk_ev,lwrk_ev,iwrk_ev,kwrk_ev,ier)
|
||||||
|
c
|
||||||
|
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,x00v,w,
|
||||||
|
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
||||||
|
. nx,tx,ny,ty,c,fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_bispev(tx,nx,ty,ny,c,kx,ky,(/xcoord0/),1,
|
||||||
|
. (/ycoord0/),1,x00,wrk_ev,lwrk_ev,iwrk_ev,kwrk_ev,ier)
|
||||||
|
c
|
||||||
|
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,y00v,w,
|
||||||
|
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
||||||
|
. nx,tx,ny,ty,c,fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_bispev(tx,nx,ty,ny,c,kx,ky,(/xcoord0/),1,
|
||||||
|
. (/ycoord0/),1,y00,wrk_ev,lwrk_ev,iwrk_ev,kwrk_ev,ier)
|
||||||
|
c
|
||||||
|
call dierckx_surfit(iopt,nisteer,xcoord,ycoord,z00v,w,
|
||||||
|
. minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps,
|
||||||
|
. nx,tx,ny,ty,c,fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier)
|
||||||
|
c
|
||||||
|
call dierckx_bispev(tx,nx,ty,ny,c,kx,ky,(/xcoord0/),1,
|
||||||
|
. (/ycoord0/),1,z00,wrk_ev,lwrk_ev,iwrk_ev,kwrk_ev,ier)
|
||||||
end if
|
end if
|
||||||
c c====================================================================================
|
|
||||||
end if
|
end if
|
||||||
c#######################################################################################
|
c#######################################################################################
|
||||||
c
|
c
|
||||||
if(fdeg.ne.0) then
|
if(incheck.ne.0) then
|
||||||
deallocate(cycoord, txycoord, cwaist1, txwaist1, cwaist2,
|
if(fdeg.ne.0) then
|
||||||
. txwaist2, crci1, txrci1, crci2, txrci2, cphi1, txphi1,
|
deallocate(c,tx,w,wrk,iwrk)
|
||||||
. cphi2, txphi2, cx0, txx0, cy0, txy0, cz0, txz0, w)
|
else
|
||||||
|
deallocate(c, tx, ty, wrk, wrk2, iwrk, w, wrk_ev,iwrk_ev,
|
||||||
|
. xpolyg, ypolyg)
|
||||||
|
end if
|
||||||
else
|
else
|
||||||
deallocate(cwaist1, txwaist1, tywaist1,
|
if(fdeg.eq.0) then
|
||||||
. cwaist2, txwaist2, tywaist2,
|
deallocate(xpolyg, ypolyg)
|
||||||
. crci1, txrci1, tyrci1, crci2, txrci2, tyrci2,
|
end if
|
||||||
. cphi1, txphi1, typhi1, cphi2, txphi2, typhi2,
|
|
||||||
. cx0, txx0, tyx0, cy0, txy0, tyy0, cz0, txz0, tyz0,
|
|
||||||
. wrk2, xpolyg, ypolyg, w)
|
|
||||||
end if
|
end if
|
||||||
c
|
c
|
||||||
c#######################################################################################
|
c#######################################################################################
|
||||||
@ -1527,8 +1486,10 @@ c set correct values for alpha, beta
|
|||||||
beta0 = ycoord0
|
beta0 = ycoord0
|
||||||
end if
|
end if
|
||||||
c#######################################################################################
|
c#######################################################################################
|
||||||
deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v,phi1v,
|
c
|
||||||
. phi2v,x00v,y00v,z00v,xcoord,ycoord)
|
c#######################################################################################
|
||||||
|
deallocate(alphastv, betastv, waist1v, waist2v, rci1v, rci2v,
|
||||||
|
. phi1v, phi2v, x00v, y00v, z00v, xcoord, ycoord)
|
||||||
c
|
c
|
||||||
return
|
return
|
||||||
end
|
end
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user