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