caf (12/09/82)
This is a C version of the MUF predicting program appearing in DEC 82
QST. It should run on v7 or usg or compatible systems. Testing so far
has consisted of comparing output with the example in the article.
with a -f flag, minimuf accepts 1 or more 10mHz readings from WWV
which are smoothed and converted to sunspot number according to the
graph in the article. Enjoy.
Chuck Forsberg, WA7KGX
Chief Engr, Computer Development Inc.
6700 S. W. 105th, Beaverton OR 97005 (503) 646-1599
...!teklabs!ogcvax!metheus!cdi!caf
/*% cc -O % -lm -o minimuf
*/
/*
* Minimuf 3.5 From Dec 1982 QST pp. 36-38.
* By Dr Paul Levine, J N Martin, D B Sailors, Robert B Rose
* Single F layer prediction Useful for 2 to 50 mHz, 250 to 6000 miles
* rms error reported to be +- 3.8 mHz.
* C version translated by Chuck Forsberg WA7KGX
*/
#include <stdio.h>
#include <math.h>
double sign();
#define PI 3.141593
double Txlat, Txlon; /* transmitter latitude, longitude */
double Rxlat, Rxlon; /* receiver latitude, longitude */
double Month; /* month, 1 to 12 */
double Day; /* day, 1 to 31 */
double Hr; /* time in UTC 0 to 24 */
double Muf; /* output MUF in mHz */
double Isn; /* sunspot number */
double Flux; /* 10mHz flux */
double r0,p1,r1,p0,k7,g1,k6,k5;
double k1,kstep,kend,a,b,p,q,c,w0,d,l0,k9,c0,t9,g7;
double t1,y2,k8,g0,m9;
double t,t4,g9,g8,u,g2,t6,u1;
char Months[] = " JanFebMarAprMayJunJulAugSepOctNovDec";
int m[] = {0,31,28,31,30,31,30,31,31,30,31,30,31};
main(argc,argv)
char *argv[];
{
int ret;
static getflux=0;
if (--argc > 0)
if (strcmp(*++argv, "-f")==0)
++getflux;
r0 = PI/180.;
p1 = 2.*PI;
r1 = 180./PI;
p0 = PI/2.;
do {
fprintf(stderr, "Transmitter Latitude Longtitude\n");
ret=scanf("%lf%lf", &Txlat, &Txlon);
if (ret==EOF)
exit(0);
}
while (ret!=2 || Txlat>90. || Txlat<-90. || Txlon>360. || Txlon<-360.);
do {
fprintf(stderr, "Receiver Latutide Longtitude\n");
ret=scanf("%lf%lf", &Rxlat, &Rxlon);
if (ret==EOF)
exit(0);
}
while (ret!=2 || Rxlat>90. || Rxlat<-90. || Rxlon>360. || Rxlon<-360.);
do {
fprintf(stderr, "Day Month : ");
ret = scanf("%lf%lf", &Day, &Month);
if (ret==EOF)
exit(0);
}
while (ret != 2 || Month<1 || Month>12 || Day<1 || Day>m[(int)Month]);
if (!getflux) {
do {
fprintf(stderr, "Sunspot number : ");
ret=scanf("%lf", &Isn);
if (ret==EOF)
exit(0);
}
while (ret != 1 || Isn <= 0);
}
else
{
fprintf(stderr,
"Daily Flux reading[s], oldest first, terminate with 0\n");
ret=scanf("%lf", &Flux);
for (;;) {
ret=scanf("%lf", &Isn);
if (ret != 1 || Isn==0)
break;
Flux = 0.8*Flux + 0.2*Isn;
}
if (Flux>130.)
Isn = (Flux-65.)*1.34;
else
Isn = Flux - 35.;
}
printf("\n\nMinimuf 3.5 Maximum Usable Frequency Prediction\n\n");
printf("Date %3.3s %2d\n", &Months[3*(int)Month], (int) Day);
printf("Transmitter Latitide %lg Longtitude %lg\n", Txlat, Txlon);
printf("Receiver Latitide %lg Longtitude %lg\n", Rxlat, Rxlon);
printf("Sunspot Number %lg\n\n", Isn);
printf(" Hour MUF(mHz)\n\n");
Txlat *= r0; Txlon *= r0; Rxlat *= r0; Rxlon *= r0;
for (Hr =0; Hr <= 23; Hr += 1.) {
minimuf();
printf(" %02d %4.1lf\n", (int)Hr, Muf);
}
exit(0);
}
minimuf()
{
k7 = sin(Txlat)*sin(Rxlat)+cos(Txlat)*cos(Rxlat)*cos(Rxlon-Txlon);
if (k7 > 1.)
k7 = 1.;
else if (k7 < -1.)
k7 = -1.;
g1 = acos(k7);
k6 = 1.59 * g1;
if (k6 < 1.)
k6 = 1.;
k5 = 1./k6;
Muf = 100.;
kstep = .9999-1./k6;
kend = 1.-1./(2.*k6);
for (k1=1./(2.*k6); kstep>0? k1 <= kend : k1 >= kend; k1 += kstep) {
if (k5 != 1.0)
k5 = 0.5;
p = sin(Rxlat);
q = cos(Rxlat);
a = (sin(Txlat)- p * cos(g1)) / (q * sin(g1));
b = g1*k1;
c = p * cos(b)+q*sin(b)*a;
d = (cos(b) - c*p) / (q*sqrt(1-(c*c)));
if (d > 1.)
d = 1.;
else if (d < -1.)
d = -1.;
d = acos(d);
w0 = Rxlon+sign(sin(Txlon-Rxlon))*d;
if (w0 < 0)
w0 += p1;
if ( w0 > p1)
w0 -= p1;
if ( c < -1.)
c = -1.;
else if (c > 1.)
c = 1.;
l0 = p0-acos(c);
t1 = 0.0172*(10.+(Month-1.)*30.4+Day);
y2 = 0.409*cos(t1);
k8 = 3.82 * w0+12.+0.13*(sin(t1)+1.2*sin(2.*t1));
k8 = k8-12.*(1.+sign(k8-24.))*sign(fabs(k8-24.));
if (cos(l0+y2) <= -0.26) {
k9 = g0 = 0;
m9 = 2.5 * g1 * k5;
if (m9 > p0)
m9 = p0;
m9 = sin(m9);
m9 = 1. + 2.5 * m9 * sqrt(m9);
}
else
{
k9 = (-0.26+sin(y2)*sin(l0))/(cos(y2)*cos(l0)+0.001);
k9 = 12. - atan(k9/sqrt(fabs(1.-k9*k9)))*7.639437;
t = k8-k9/2.+Rxlat*(1.-sign(k8-k9/2.))*sign(fabs(k8-k9/2.));
t4 = k8+k9/2.-12.*(1.+sign(k8+k9/2.-24.))*sign(fabs(k8-k9/2.-24.));
c0 = fabs(cos(l0+y2));
t9 = 9.7*pow(c0,9.6);
if (t9 < 0.1)
t9 = 0.1;
m9 = 2.5*g1*k5;
if (m9 > p0)
m9 = p0;
m9 = sin(m9);
m9 = 1. + 2.5*m9*sqrt(m9);
if ((t4 < t && (Hr-t4)*(t-Hr) > 0)
|| (t4 >= t && (Hr-t)*(t4-Hr) <= 0))
{
t6 = Hr+12.*(1.+sign(t4-Hr))*sign(fabs(t4-Hr));
g8 = PI*t9/k9;
u = (t4-t6)/2.;
u1 = -k9/t9;
g0 = c0*(g8*(exp(u1)+1.))*exp(u)/(1.+g8*g8);
}
else
{
t6 = Hr+12.*(1.+sign(t-Hr))*sign(fabs(t-Hr));
g9 = PI * (t6-t)/k9;
g8 = PI * t9/k9;
u = (t-t6)/t9;
g0 = c0*(sin(g9)+g8*(exp(u)-cos(g9)))/(1.+g8*g8);
g7=c0*(g8*(exp(-k9/t9)+1.))*exp((k9-24.)/2.)/(1.+g8*g8);
if (g0 < g7)
g0 = g7;
}
}
g2 =(1.+Isn/250.)*m9*sqrt(6.+58.*sqrt(g0));
g2 = g2*(1.-0.1*exp((k9-24.)/3.));
g2 = g2*(1.+(1.-sign(Txlat)*sign(Rxlat))*0.1);
g2 = g2*(1.-0.1*(1.+sign(fabs(sin(l0))-cos(l0))));
if (Muf > g2)
Muf = g2;
}
return;
}
double
sign(g)
double(g);
{
if (g > 0.)
return 1.;
if (g < 0.)
return -1.;
return 0.0;
}