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