subroutine ray_qp
c
c s/r to find the values of X and its partial derivatives pxpr,
c pxpth and pxpr to be used with the Jones program. The parameters
c of the ionospheric model have previously been set up at a special
c grid by the program ION_QP. This grid has a north pole at the TX.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c it is assumed that raytracing is to be done with the earth's
c magnetic field included --- w(1) = +/- 1.
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c written by leo mcnamara 31-jul-90
c modified 3-jun-91 to allow for cases with grid longitudes around 0
c----- last modified 3-jun-91
c
c USES s/rs QP_XPXPR, XPXPR_QPAR, QLFNH
c
logical lon_flag
character*8 modx,id
real *8 x,freq,hh,pxpr,pxpth,pxpph,pxpt,hmax,r,w0,w
c
common / xx / modx(2),x,pxpr,pxpth,pxpph,pxpt,hmax
common r(6) / ww / id(10),w0,w(400)
c
dimension p(50,5,20),kr(4),kc(4),x_colat(4),y_lon(4),xs(4),
a pr(4),hm(4),colat(50),glon(5)
c
data pi,radeg,degrad / 3.141592654,57.29577951,0.017453292 /
data ipass / 0 /
c
c-----------------------------
entry ELECTX
c-----------------------------
c
if(w(1).eq.0.) stop ' Must have magnetic field '
c
c read in the ionospheric model on the first pass
if(ipass.ne.0) go to 100
lon_flag = .false.
open (unit=3,file='ionmod_qp.dat',type='old')
c read in 5 header lines
do 1 i = 1,5
1 read(3,333) head
333 format(a4)
c read in the coords of the TX after it has been moved 1 km down range
read(3,102) tlat_moved,tlon_moved
c overwrite any specified TX location
w(4) = dble(tlat_moved) * degrad
w(5) = dble(tlon_moved) * degrad
c read in the geom coords of the TX
read(3,102) tx_gmlat,tx_gmlon
102 format(5e16.8)
c the North pole is at the geomag coords of the TX
pole_lat = tx_gmlat
pole_lon = tx_gmlon
c read in the co-latitude values of the latitude circles of the grid
read(3,101) nlat
101 format(i4)
read(3,102) (colat(i),i=1,nlat)
do 1101 i = 1,nlat
1101 print*,' colat ', i,colat(i)
c read in the longitudes of the grid
read(3,101) nlon
read(3,102) (glon(i),i=1,nlon)
c See if the longitudes straddle zero degrees. If so, set lon_flag
c to TRUE and add 360 to the easterly longitudes.
if(glon(1).gt.340.0.and.glon(nlon).lt.20.) lon_flag = .true.
do 1102 i = 1,nlon
if(lon_flag.and.glon(i).lt.20.0) glon(i) = glon(i) + 360.
1102 print*,' glon ', i,glon(i)
c
c search relies on the COLAT and GLON increasing monotonically - check
c that this is true. Could have problems if longitude of RX in GEOTRANS
c coords is around 0 (in spite of fix above! - let's be careful)
do 104 i = 2,nlat
if(colat(i).le.colat(i-1)) stop ' Problems with COLAT grid'
104 continue
do 103 i = 2,nlon
print*,' glon(i),glon(i-1) ', i,glon(i),glon(i-1)
if(glon(i).le.glon(i-1)) stop ' Problems with LON grid'
103 continue
c
c read in the parameters of the ionospheric model at each gridpoint
c there are NPARAM parameters at each grid point
read(3,101) nparam
print*,' NPARAM,nlat,nlon ', nparam,nlat,nlon
do 105 i = 1,nlat
do 105 j = 1,nlon
105 read(3,102) (p(i,j,l),l = 1,nparam)
c
c geographic coords of TX
gglat_tx = p(1,1,1)
gglon_tx = p(1,1,2)
print*,' Geog coords of TX ',gglat_tx,gglon_tx
c
ipass = 1
c return on first call, to avoid problems with rest of s/r
return
c
100 continue
c
freq = w(6)
hh = r(1) - w(2)
c coords of point - r,theta,phi - work in degrees
rr = sngl(r(1))
theta = sngl(r(2)) * radeg
phi = sngl(r(3)) * radeg
if(phi.lt.0.) phi = phi + 360.
gmlat = 90. - theta
gmlon = phi
ccc print*,' RAY_QP rr, theta,phi ', rr-6370.,theta,phi
c
c geog coords - we need this, since we cannot go directly from
c GEOM coords to GEOTRANS coords (too dumb)
gmpole_lat = 78.6
gmpole_lon = 290.2
call gg_gm (-1,gmpole_lat,gmpole_lon,gglat,gglon,90.-theta,phi)
ccc print*,' geog coords ',rr,gglat,gglon
c
c return if point is clearly below the ionosphere
x = 0.
pxpr = 0.
pxpth = 0.
pxpph = 0.
pxpt = 0.
if(hh.lt.50.) return
c
c coords of point in GEOTRANS coord frame (North pole at TX)
c work from geog coords (since we can't work out how to go directly
c from geom to geotrans coords)
call gg_gm (1,gglat_tx,gglon_tx,gglat,gglon,xlat,xlon)
c modify xlon if the longitudes straddle zero
if(lon_flag.and.xlon.lt.20.0) xlon = xlon + 360.
ccc write(6,666) hh,90.-xlat,xlon
666 format(2x,'GEOTRANS location ', 3f10.2)
c
c find the gridpoints surrounding (xlat,xlon)
c co-latitude of the point (xlat,xlon)
xcolat = 90. - xlat
nlatm1 = nlat - 1
do 200 i=1,nlatm1
iclat = i
if(xcolat.gt.colat(i).and.xcolat.le.colat(i+1)) go to 222
200 continue
print*,' xcolat,colat(1),colat(nlat) ',
a xcolat,colat(1),colat(nlat)
stop ' Point outside grid values of co-latitude'
c
222 continue
nlonm1 = nlon - 1
do 250 i = 1,nlonm1
ilon = i
if(xlon.gt.glon(i).and.xlon.le.glon(i+1)) go to 226
250 continue
print*,' xlon,glon(1),glon(nlon) ',
a xlon,glon(1),glon(nlon)
stop ' Point outside grid values of longitude'
226 continue
c
ccc print*,' iclat,ilon ', iclat,ilon
c
c point is surrounded by the 4 gridpoints defined by colat(iclat)
c and colat(iclat+1) and glon(ilon) and glon(ilon+1)
kr(1) = iclat
kr(2) = iclat
kr(3) = iclat + 1
kr(4) = iclat + 1
kc(1) = ilon
kc(2) = ilon + 1
kc(3) = ilon + 1
kc(4) = ilon
c
c Interpolate to get values of X, pX/pr and hmF2
do 300 k = 1,4
c set up model parameters at each grid point k1,k2 (row,col)
k1 = kr(k)
k2 = kc(k)
foe = p(k1,k2,5)
yme = p(k1,k2,6)
hme = p(k1,k2,7)
fof1 = p(k1,k2,8)
ymf1 = p(k1,k2,9)
hmf1 = p(k1,k2,10)
fof2 = p(k1,k2,11)
ymf2 = p(k1,k2,12)
hmf2 = p(k1,k2,13)
hsf = p(k1,k2,14)
hsf2 = p(k1,k2,15)
fnmax = p(k1,k2,16)
ymax = p(k1,k2,17)
himax = p(k1,k2,18)
invqp = p(k1,k2,19)
c
c calculate the values of X and pX/pr at each gridpoint
call qp_xpxpr
a (foe,yme,hme,fof1,ymf1,hmf1,fof2,ymf2,hmf2,hsf,hsf2,
a fnmax,ymax,himax,invqp)
c x, pxpr and hmf2 are in labelled common /XX/
c
xs(k) = sngl(x)
pr(k) = sngl(pxpr)
hm(k) = hmf2
300 continue
c
ccc print*,' x ', (xs(k),k=1,4)
c return if all 4 values of X = 0
if(xs(1)+xs(2)+xs(3)+xs(4).lt.1.e-10) return
c
c interpolate at the point P to get x, hmax and pxpr (2-D interpolation)
c the lat and lon of the 4 corners
x_colat(1) = colat(iclat)
x_colat(2) = colat(iclat)
x_colat(3) = colat(iclat + 1)
x_colat(4) = colat(iclat + 1)
y_lon(1) = glon(ilon)
y_lon(2) = glon(ilon + 1)
y_lon(3) = glon(ilon + 1)
y_lon(4) = glon(ilon)
ccc print*,' XC ',(x_colat(i),i=1,4)
ccc print*,' YC ',(y_lon(i),i=1,4)
call interp_2d (xcolat,xlon,x_colat,y_lon,xs,ans_x)
call interp_2d (xcolat,xlon,x_colat,y_lon,pr,ans_px)
call interp_2d (xcolat,xlon,x_colat,y_lon,hm,ans_hm)
x = dble(ans_x)
ccc print*,' Interpolated X ', x
pxpr = dble(ans_px)
hmax = dble(ans_hm)
c
c calculate the gradients wrt theta and phi, by finding X at a nearby
c point P1 with (1) theta + DELTA_theta, same phi, ==> pxpth
c point P2 with (2) theta, phi + DELTA_phi ==> pxpph
c The steps DELTA are made small enough that it is REASONABLY safe
c to use the same interpolation grid as before
delta = 0.2
c we need the size of the interpolation box in geom coords
size_gmlat = abs(p(iclat,ilon,3) - p(iclat+1,ilon,3))
c all the longitudes are identical at the TX (iclat = 1)
jclat = iclat
if(iclat.eq.1) jclat = 2
size_gmlon = abs(p(jclat,ilon,4) - p(jclat,ilon+1,4))
delta_theta = delta * size_gmlat
delta_phi = delta * size_gmlon
c find GEOTRANS coords of the two points, from the geom coords,
c via the geog coords
gmlat1 = gmlat + delta_theta
gmlon1 = gmlon
gmlat2 = gmlat
gmlon2 = gmlon + delta_phi
call gg_gm
a (-1,gmpole_lat,gmpole_lon,gglat1,gglon1,gmlat1,gmlon1)
call gg_gm
a (-1,gmpole_lat,gmpole_lon,gglat2,gglon2,gmlat2,gmlon2)
call gg_gm (1,gglat_tx,gglon_tx,gglat1,gglon1,xlat1,xlon1)
call gg_gm (1,gglat_tx,gglon_tx,gglat2,gglon2,xlat2,xlon2)
c
xcolat1 = 90. - xlat1
xcolat2 = 90. - xlat2
ccc print*,' GEOTRANS 1 & 2 ', xcolat1,xlon1,xcolat2,xlon2
ccc print*,' x ', (xs(k),k=1,4)
call interp_2d (xcolat1,xlon1,x_colat,y_lon,xs,ans_x1)
call interp_2d (xcolat2,xlon2,x_colat,y_lon,xs,ans_x2)
c
pxpth = dble(ans_x1 - ans_x) / (delta_theta * degrad)
c the gradient wrt co-latitude = - gradient wrt latitude
pxpth = - pxpth
pxpph = dble(ans_x2 - ans_x) / (delta_phi * degrad)
c
ccc print*, ans_x,ans_x1,ans_x2,pxpth,pxpph
c
c there is no variation with time
pxpt =0.
c
return
end
subroutine interp_2d (xp,yp,xc,yc,z,zp)
dimension xc(1),yc(1),z(1)
xa = xc(1)
xd = xc(4)
ya = yc(1)
yb = yc(2)
if(xa.eq.xd) type*,'theta,phi ',th,ph
if(ya.eq.yb) type*,'theta,phi ',th,ph
if(xa.eq.xd) type*,' xa & xd are equal --- ',xa
if(ya.eq.yb) type*, 'ya & yb are equal --- ',ya
if(xa.eq.xd) type*,' xc', (xc(i),i=1,4)
if(xa.eq.xd) type*,' yc', (yc(i),i=1,4)
if(ya.eq.yb) type*,' xc ',(xc(i),i=1,4)
if(ya.eq.yb) type*,' yc ',(yc(i),i=1,4)
za = z(1)
zb = z(2)
zc = z(3)
zd = z(4)
fac1 = (yp - ya) / (yb - ya)
fac2 = (xp - xa) / (xd - xa)
zp = za + fac1*(zb-za) + fac2*(zd-za+fac1*(za-zb+zc-zd))
return
end