c..... file name dudtilt.for
SUBROUTINE DUDTILT(ipass,ff)
C DUDENEY N(H) PROFILE S/R FOR USE WITH Jones RAYTRACING PROG
C THAT CONSIDERS A TILTED (july 87)IONOSPHERE
c calculate the parameters x,pxpr,pxpth,pxpph corresponding to a
c given position above the surface of the earth, for a wave freq ff.
C SEE JATP 40(2), 1978, PG 195
C LEO MCNAMARA APRIL 1981
c----- LAST MODIFIED 28-sep-87
c
C.....REQUIRES S/R DIMSOL & FUNCTIONS HMAX,YMAX (DUDENEY S/RS)
C
C S/R IS ESSENTIALLY THE SAME AS DUDENEY'S FUNCTION FN(H), EXCEPT
C THAT THE SLOPE PXPR IS ALSO CALCULATED
C........................................................................
C
c IPASS = 0 to calculate Dudeney parameters; set=1 afterwards
C........................................................................
C FOF2 F2 LAYER CRITICAL MANDATORY
C FOF1 F1 LAYER CRITICAL NO F1 GENERATED IF <=0
C FOE E LAYER CRITICAL SET=.6 IF <=0
C HFF2 MIN F2 VIRTUAL HEIGHT MANDATORY IF YMF2<=0
C M3000 M(3000)F2 MANDATORY IF HMF2<=0
C HME E REGION PEAK HEIGHT SET=100 IF <=0
C YME E REGION SEMITHICKNESS SET=20 IF <=0
C HMF2 F REGION PEAK HEIGHT CALCULATED IF <=0
C YMF2 F REGION SEMITHICKNESS CALCULATED IF <=0
C H0 BASE HEIGHT OF IONOSPHERE E LAYER IF XD<=0
C XD D REGION TANGENT TO E LAYER AT XD*FOE
C DV E-F FRACTIONAL VALLEY DEPTH
C........................................................................
C
REAL M3000,h
real*8 x,pxpr,pxpth,pxpph,pxpt,hmax,ff,hh,r,w0,w
character*8 modx,id
C
C DUDENEY'S LABELLED COMMON - OUTPUT FROM DIMSOL
COMMON / PARMS / BETA,HT,FT,A,HF,AD,TAU,HD
C JONES CODE
COMMON / XX / MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX
common r(6) /ww/ id(10),w0,w(400)
c
common / iondat / fof2,fof1,foe,hff2,m3000,hme,yme,
a hmf2,ymf2,h0,xd,dv
c
common / ioncent / hapog, tilt, cosba, sinba,
& cosblon, sinblon, cosblat, sinblat
c
data modx / 'DUDNEY',' ' /
DATA PI,C1 / 3.1415926536,1. /
C
fn=0.
X=0.d0
PXPR=0.d0
PXPTH=0.d0
PXPPH=0.d0
f = sngl (ff)
c
C FIRST CALL TO S/R - IPASS=0 - CALCULATE DUDENEY PROFILE
IF(IPASS.EQ.1) GO TO 10
CALL DIMSOL(IERROR)
HMAX= dble(HMF2)
IF(IERROR.EQ.0) GO TO 10
PRINT 200,IERROR
200 FORMAT(//5X,'*DIMSOL RETURNED ERROR CODE*',I3//)
STOP 200
10 continue
c
c-------------------------------------------------------------------------------
c
c in this section, created 3-jul-87, the profile above is considered
c to be a function only of distance about some point (xion,yion,zion).
c This point depends on the tilt we assign to it which is derived
c from a measured "tilt" by the subroutine "lr_tilt".
c The distance from this point to the beacon that measures tilt
c (see beacontilt) is w(2), the radius of the earth.
c given r(1), r(2), r(3) (=r,co-lat,long) of the point we consider
c we calculate a new radius, r_dash.
c Later we evaluate X and pxpr_dash and
c from co-ordinate transforms and the chain rule evaluate
c pxpr, pxpth, pxpph. (david barker)
c
c calculate the cartesian co-ords of the point r
r1 = sngl (r(1))
r2 = sngl (r(2))
r3 = sngl (r(3))
ccc type*,' r1,r2,r3 ',r1,r2,r3
if(r3.lt.0.0) r3 = r3 + 2.*pi
cosr2 = cos(r2)
cosr3 = cos(r3)
sinr2 = sin(r2)
sinr3 = sin(r3)
xi = r1 * sinr2 * cosr3
yi = r1 * sinr2 * sinr3
zi = r1 * cosr2
c
c calculate the applicable ionospheric tilt "tilth"
hh = r(1) - w(2)
h = sngl(hh)
call lr_tilt(tilt,hapog,h,tilth)
c
c calculate the iono-centre applicable to the tilt "tilth"
sth = sin(tilth)
cth = cos(tilth)
w2 = sngl(w(2))
xion = w2 * ( cosblon*cosblat*(1.-cth)
& +sth*(sinblon*sinba
& +cosblon*sinblat*cosba))
yion = w2 * ( sinblon*cosblat*(1.-cth)
& -sth*(cosblon*sinba
& -sinblon*sinblat*cosba))
zion = w2 * ( sinblat*(1.-cth)
& -cosblat*sth*cosba)
c
c calculate the new r co-ord in the new reference frame
r_dash = sqrt((xi-xion)*(xi-xion) +
& (yi-yion)*(yi-yion) +(zi-zion)*(zi-zion))
c
hh = dble(r_dash) -w(2)
h = sngl(hh)
ccc type*,' w(2),r_dash,hh ', w(2),r_dash,hh
c
c-------------------------------------------------------------------------------
C
C CHECK TO SEE IF POINT IS BELOW BASE HEIGHT
IF(H.LE.H0) RETURN
C
C D-REGION SEGMENT
IF(H.GT.HD) GO TO 20
FN=TAU/ALOG(AD/(H-H0))
DFNDH=FN*FN/TAU/(H-H0)
GO TO 60
C
C E-REGION SEGMENT DUDENEY EQN 5
20 IF(H.GT.HME) GO TO 25
CON1=(HME-H)/YME
CON2=SQRT(C1-CON1*CON1)
FN=FOE*CON2
DFNDH=FOE*CON1/CON2/YME
GO TO 60
c
25 continue
C
C F1-REGION SEGMENT APPROX DUDENEY EQN 12
IF(H.GT.HT) GO TO 40
ANG1=2.*PI*(H-HME)/(HT-HME)
ANG2=(H-HME)/BETA/YMF2
COS1=COS(ANG1)
SIN1=SIN(ANG1)
COS2=COS(ANG2)
SIN2=SIN(ANG2)
FN=FOE*(A*(C1-COS1)+C1/COS2)
DFNDH=2.*PI*A/(HT-HME)*SIN1+C1/BETA/YMF2/COS2/COS2*SIN2
DFNDH=FOE*DFNDH
GO TO 60
C
C F2-REGION DUDENEY EQN 1
C TOPSIDE DENSITY DROPS TO ZERO AT H=HMF2+(PI/2)*YMF2
40 IF(H.GT.HMF2+(PI/2.)*YMF2) RETURN
ANG=(HMF2-H)/YMF2
FN=FOF2*COS(ANG)
DFNDH=FOF2*SIN(ANG)/YMF2
C
C CONVERT FROM FN TO X AND FROM DFNDH TO PXPR (for Jones prog)
60 F2=F*F
X= dble (FN*FN/F2)
PXPR= dble (2.*FN*DFNDH/F2)
C
c-------------------------------------------------------------------------------
c
c in this section we calculate pxpr proper (the above code calculated
c pxpr_dash) and the gradients pxpth, pxpph that we created by tilting
c the ionosphere
c
c calculate conversion factors
convr = (r1 -xion * sinr2 * cosr3
& -yion * sinr2 * sinr3
& -zion * cosr2 ) / r_dash
convth = r1 * (-xion * cosr2 * cosr3
& -yion * cosr2 * sinr3
& +zion * sinr2 ) / r_dash
convph = r1 * ( xion * sinr2 * sinr3
& -yion * sinr2 * cosr3 ) / r_dash
c
ccc type*,' pxpr,convr,convth,convph=',pxpr,convr,convth,convph
c now calculate the derivatives
pxpph = pxpr * dble(convph)
pxpth = pxpr * dble(convth)
pxpr = pxpr * dble(convr)
ccc type*, ' pxpr,pxpth,pxpph=', pxpr,pxpth,pxpph
c
c-------------------------------------------------------------------------------
c
RETURN
END