/* *********************************************************************** * Program to calculate maximum usable frequency and received signal * * strength for high-frequency radio circuits. Uses MINIMUF 3.5. * *********************************************************************** input format: first line contains five numbers: 3 1 1 194 20 1. output format 1 receiver power (dB/uV) 2 elevation angle (deg) 3 delay (ms) 2. month (1-12) 3. day (1-31) 4. smoothed sunspot number 5. transmitter ERP (dBW) second line contains transmitter coordinates and name: 39.7000 -75.7825 W3HCF Newark third and subsequent lines contain receiver coordinates and name: 39.6808 -75.7486 Evans Hall 45.18 -75.45 CHU Ottawa 40.6803 -105.0408 WWV Ft Collins 21.9906 -159.7667 WWVH Hawaii 52.3667 -1.1833 MSF Rugby timecode transmitters, frequencies and geographic coordinates WWV Ft. Collins, CO 2.5, 5, 10, 15, 20 MHz 40:40N, 105:03W WWVB Ft. Collins, CO 60 kHz 40:40:49N, 105:02:27W WWVH Kauai, HI 2.5, 5, 10, 15, 20 MHz 21:59:26N, 159:46:00W CHU Ottawa, Canada 3330, 7335, 14670 45:18N, 75:45N MSF Rugby, UK 60 kHz, 2.5, 5, 10 MHz 52:22N, 1:11W */ #include #include #include /* Parameters */ #define R 6371.2 /* radius of the Earth (km) */ #define hE 110. /* height of E layer (km) */ #define hF 320. /* height of F layer (km) */ #define GAMMA 1.42 /* geomagnetic constant */ #define PI 3.141592653589 /* the real thing */ #define PIH PI/2. /* the real thing/2 */ #define PID 2.*PI /* the real thing*2 */ #define VOFL 2.9979250e8 /* velocity of light */ #define D2R PI/180. /* degrees to radians */ #define R2D 180./PI /* radians to degrees */ #define MINBETA 5.*D2R /* min elevation angle (rad) */ #define RINZ 50. /* receiver input impedance (ohms) */ #define BOLTZ 1.380622e-23 /* Boltzmann's constant */ #define NTEMP 273. /* receiver noise temperature (deg K) */ #define DELTAF 2400. /* communication bandwidth (Hz) */ #define MPATH 3. /* multipath threshold (dB) */ #define GLOSS 3. /* ground-reflection loss (dB) */ #define FMAX 8 /* max frequencies */ #define HMAX 10 /* max hops */ /* Function declarations */ extern FILE *fopen(); static double minimuf(void), sgn(double); void edit(int); /* Global declarations */ int month; /* month of year (1-12) */ int day; /* day of month*/ int hour; /* hour of day (UTC) */ double freq[FMAX]; /* frequencies (MHz) */ double gain[46][FMAX]; /* antenna gain (main lobe) (dB) */ double dB1; /* transmitter output power (dBW) */ double ssn; /* smoothed sunspot number */ double lat1; /* transmitter latitude (deg N) */ double lon1; /* transmitter longitude (deg W) */ char site1[30]; /* transmitter site name */ double lat2; /* receiver latitude (deg N) */ double lon2; /* receiver longitude (deg W) */ char site2[30]; /* receiver site name */ double b1; /* transmitter bearing (rad) */ double b2; /* receiver bearing (rad) */ double d; /* great-circle distance (rad) */ double theta; /* hour angle (rad) */ int hop; /* number of ray hops */ double phi; /* angle of incidence (rad) */ FILE *fp_in; /* file handle */ char fn_in[12] = "gain.dat"; /* antenna file name */ char fq_in[12] = "ntp.dat"; /* qth file name */ int flag; /* output format */ /* Main program */ main() { /* Hour variables */ int offset; /* offset for local time (hours) */ int time; /* local time at receiver */ char daynight[HMAX]; /* day/night indicator (jnx) */ double mufE[HMAX]; /* max E-layer muf (MHz) */ double mufF[HMAX]; /* min F-layer muf (MHz) */ double absorp[HMAX]; /* ionospheric absorption coefficient */ double dB2[HMAX]; /* receiver input power (dBuV) */ double path[HMAX]; /* path length (km) */ double beta[HMAX]; /* elevation angle (rad) */ double ant[HMAX]; /* antenna gain (dB) */ /* Path variables */ double delay; /* path distance (ms) */ double dhop; /* hop angle/2 (rad) */ int daytime, nightime; /* path indicators */ double psi; /* sun zenith angle (rad) */ double xr, yr; /* reflection zone coordinates (rad) */ double xs, ys, yp; /* subsolar coordinates (rad) */ double fcE; /* E-layer critical frequency (MHz) */ double fcF; /* F-layer critical frequency (MHz) */ char cutoff; /* layer cutoff indicator (EFM) */ double beta1, phi1, level, muf, dist, Z, floor; /* double temporaries */ int i,j,h,n; /* int temporaries */ float fp1, fp2; /* float temporaries */ /* Establish initial conditions */ if ((fp_in = fopen (fn_in, "r")) == NULL) { printf ("Antenna file %s not found\n", fn_in); exit (1); } for (j = 0; j < FMAX; j++) { fscanf(fp_in, " %f ", &fp1); freq[j] = fp1; } for (i = 0; i < 46; i++) { for (j = 0; j < FMAX; j++) { fscanf(fp_in, " %f ", &fp1); gain[i][j] = fp1; } } if ((fp_in = fopen (fq_in, "r")) == NULL) { printf ("QTH file %s not found\n", fq_in); exit (2); } if (fscanf(fp_in, " %i %i %i %f %f", &flag, &month, &day, &fp1, &fp2) != 5) exit (0); ssn = fp1; dB1 = fp2; if (fscanf(fp_in, " %f %f %[^\n]", &fp1, &fp2, site1) != 3) exit (0); lat1 = fp1*D2R; lon1 = -fp2*D2R; printf("Transmitter %s\n", site1); L1: if (fscanf(fp_in, " %f %f %[^\n]", &fp1, &fp2, site2) != 3) exit (0); lat2 = fp1*D2R; lon2 = -fp2*D2R; printf("Receiver %s\n", site2); /* Compute great-circle bearings, great-circle distance, min hops and path delay */ theta = lon1-lon2; if (theta >= PI) theta = theta-PID; if (theta <= -PI) theta = theta+PID; d = acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(theta)); if (d < 0.) d = PI+d; b1 = acos((sin(lat2)-sin(lat1)*cos(d))/(cos(lat1)*sin(d))); if (b1 < 0.) b1 = PI+b1; if (theta < 0) b1 = PID-b1; b2 = acos((sin(lat1)-sin(lat2)*cos(d))/(cos(lat2)*sin(d))); if (b2 < 0.) b2 = PI+b2; if (theta >= 0.) b2 = PID-b2; hop = d/(2.*acos(R/(R+hF))); beta1 = 0.; while (beta1 < MINBETA) { hop = hop+1; dhop = d/(hop*2.); beta1 = atan((cos(dhop)-R/(R+hF))/sin(dhop)); } phi = R*cos(beta1)/(R+hE); phi = atan(phi/sqrt(1.-phi*phi)); delay = 2.*hop*sin(dhop)*(R+hF)/cos(beta1)/VOFL*1e6; printf("\nSmoothed sunspot number:%4.0f Month:%3i Day:%3i\n", ssn, month, day); printf("Power:%3.0f dBW Distance:%6.0f km Delay:%5.1f ms\n", dB1, d*R, delay); printf("Location Lat Long Azim\n"); printf("%-27s %7.2fN %7.2fW %3.0f\n", site1, lat1*R2D, lon1*R2D, b1*R2D); printf("%-27s %7.2fN %7.2fW %3.0f\n", site2, lat2*R2D, lon2*R2D, b2*R2D); printf("UT LT MUF%8.1f%8.1f%8.1f%8.1f%8.1f%8.1f%8.1f%8.1f\n" , freq[0], freq[1], freq[2], freq[3], freq[4], freq[5], freq[6], freq[7]); /* Hour loop: This loop determines the min-hop path and next two higher-hop paths. It selects the most likely path for each frequency and calculates the receiver power. The minimum F-layer MUF is computed directly from MINIMUF 3.5 and regardless of the actual number of hops, geographic coordinates or time of day. */ offset = lon2*12./PI+.5; for (hour = 0; hour < 24; hour++) { time = hour-offset; if (time < 0) time = time+24; if (time >= 24) time = time-24; muf = minimuf(); fcF = muf * cos(phi); printf("%2i %2i", hour, time); /* subsolar coordinates */ xs = (month-1.)*365.25/12.+day-80.; xs = 23.5*D2R*sin(xs/365.25*PID); ys = (hour*15.-180.)*D2R; /* Path loop: This loop determines the geometry of the min-hop path and the next two hiher-hop paths. It calculates the minimum F-layer MUF, maximum E-layer MUF and ionospheric absorption factor for each path. */ for (h = hop; h < hop+2; h++) { daytime = 0.; nightime = 0.; mufE[h] = 0.; absorp[h] = 0.; daynight[h] = ' '; /* path geometry, min F-layer MUF */ dhop = d/(h*2.); beta1 = atan((cos(dhop)-R/(R+hF))/sin(dhop)); phi1 = R*cos(beta1)/(R+hE); phi1 = atan(phi1/sqrt(1.-phi1*phi1)); mufF[h] = fcF/cos(phi1); /* Hop loop: This loop calculates the reflection zones and subsolar zenith angle for each hop along the path. It then computes the minimum E-layer MUF and ionospheric absorption coeeficient for the total path. */ for (dist = dhop; dist < d; dist = dist+dhop*2) { /* reflection zone coordinates */ xr = acos(cos(dist)*sin(lat1)+sin(dist)*cos(lat1)*cos(b1)); if (xr < 0.) xr = PI+xr; xr = PIH-xr; yr = acos((cos(dist)-sin(xr)*sin(lat1))/(cos(xr)*cos(lat1))); if (yr < 0.) yr = PI+yr; if (theta < 0.) yr = -yr; yr = lon1-yr; if (yr >= PI) yr = yr-PID; if (yr <= -PI) yr = yr+PID; yp = ys-yr; if (yp > PI) yp = PID-yp; if (yp < -PI) yp = -PID+yp; /* E-layer MUF */ psi = acos(sin(xr)*sin(xs)+cos(xr)*cos(xs)*cos(yp)); if (psi < 0.) psi = PI+psi; Z = cos(psi); if (Z > 0.) { Z = .9*pow((180.+1.44*ssn)*Z,.25); if (Z < .005*ssn) Z = .005*ssn; } else Z = .005*ssn; Z = Z/cos(phi1); if (Z > mufE[h]) mufE[h] = Z; /* ionospheric absorption coeeficient */ Z = psi; if (Z > 100.8*D2R) { Z = 100.8*D2R; nightime++; } else daytime++; Z = cos(90./100.8*Z); if (Z < 0.) Z = 0; Z = (1.+.0037*ssn)*pow(Z,1.3); if (Z < .1) Z = .1; absorp[h] = absorp[h]+Z; } /* path flags */ if (daytime > 0.) daynight[h] = 'j'; if (nightime > 0.) daynight[h] = 'n'; if ((daytime > 0.) && (nightime > 0.)) daynight[h] = 'x'; } /* Frequency loop: This loop calculates the receiver power for each frequency and path. It discards paths above the minimum F-layer MUF or below the maximum E-layer MUF and selects the path with maximum power. It then calculates the combined power of the remaining paths and determines if multipath conditions exist. */ printf("%5.1f ", mufF[hop]); for (i = 0; i < FMAX; i++) { level = -1e6; j = 0; /* receiver power for each path */ for (h = hop; h < hop+2; h++) { dhop = d/(h*2.); switch (daynight[h]) { case 'n': Z = 250.; break; case 'j': Z = 350.; break; default: Z = hF; } beta[h] = atan((cos(dhop)-R/(R+Z))/sin(dhop)); phi1 = R*cos(beta[h])/(R+hE); phi1 = atan(phi1/sqrt(1.-phi1*phi1)); path[h] = 2.*h*sin(dhop)*(R+Z)/cos(beta[h]); dB2[h] = -1e6; if ((freq[i] < mufF[h]) && (freq[i] > mufE[h])) { Z = dB1+137.; Z = Z-32.44-20.*log10(path[h])-h*GLOSS; Z = Z-8.686*log(freq[i]); Z = Z-677.2*absorp[h]/cos(phi1) /(pow((freq[i]+GAMMA), 1.98)+10.2); muf = modf(beta[h]*R2D/2., &dist); n = dist; ant[h] = gain[n][i]+muf*(gain[n+1][i]-gain[n][i]); Z = Z+ant[h]; dB2[h] = Z; if (dB2[h] > level) { level = dB2[h]; j = h;; } } } /* select path and determine multipath */ floor = 20.*log10(sqrt(4.*RINZ*BOLTZ*NTEMP*DELTAF))+120.; muf = 0.; if (j > 0) { for (h = hop; h < hop+2; h++) muf = muf+exp(dB2[h]/10.); muf = 10.*log10(muf); cutoff = ' '; if (dB2[j] < muf+MPATH) cutoff = 'm'; if (dB2[j] < floor) cutoff = 's'; switch (flag) { case 1: printf("%5.0f%c%1i%c", dB2[j], daynight[j], j, cutoff); break; case 2: printf("%5.1f%c%1i%c", beta[j]*R2D, daynight[j], j, cutoff); break; case 3: printf("%5.1f%c%1i%c", path[j]/300, daynight[j], j, cutoff); } } else { printf(" "); } } printf("\n"); } goto L1; } /* MINIMUF 3.5 (From QST December 1982) Inputs lat1 = transmitter latitude, lon1 = transmitter longitude lat2 = receiver latitude, lon2 = receiver longitude month = month of year day = day of month hour = UTC hour, ssn = sunspot number Outputs mufF = F-layer MUF */ double minimuf() { double MUF; double A, B, C, D, P, Q; double Y1, Y2; double T, T4, T6, T9; double G0, G1, G2, G7, G8, G9; double K, K1, K5, K6, K7, K8, K9; double M9, C0, U, U1, W0, L0; K7 = sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon2-lon1); if (K7 < -1.) K7 = -1.; if (K7 > 1.) K7 = 1.; G1 = acos(K7); K6 = 1.59*G1; if (K6 < 1.) K6 = 1.; K5 = 1./K6; MUF = 100.; for (K1 = 1./(2.*K6); K1 <= 1.-1./(2.*K6); K1 = K1+fabs(.9999-1./K6)) { if (K5 != 1.) K5 = .5; P = sin(lat2); Q = cos(lat2); A = (sin(lat1)-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.; if (D > 1.) D = 1.; D = acos(D); W0 = lon2+sgn(sin(lon1-lon2))*D; if (W0 < 0) W0 = W0+PID; if (W0 >= PID) W0 = W0-PID; if (C < -1.) C = -1.; if (C > 1.) C = 1.; L0 = PIH-acos(C); Y1 = .0172*(10.+(month-1.)*30.4+day); Y2 = .409*cos(Y1); K8 = 3.82*W0+12.+.13*(sin(Y1)+1.2*sin(2.*Y1)); K8 = K8-12.*(1.+sgn(K8-24.))*sgn(fabs(K8-24.)); if (cos(L0+Y2) > -.26) goto s1520; K9 = 0.; G0 = 0.; M9 = 2.5*G1*K5; if (M9 > PIH) M9 = PIH; M9 = sin(M9); M9 = 1.+2.5*M9*sqrt(M9); goto s1770; s1520: K9 = (-.26+sin(Y2)*sin(L0))/(cos(Y2)*cos(L0)+.001); K9 = 12.-atan(K9/sqrt(fabs(1.-K9*K9)))*7.639437; T = K8-K9/2.+12.*(1.-sgn(K8-K9/2.))*sgn(fabs(K8-K9/2.)); T4 = K8+K9/2.-12.*(1.+sgn(K8+K9/2.-24.))*sgn(fabs(K8+K9/2-24.)); C0 = fabs(cos(L0+Y2)); T9 = 9.7*pow(C0,9.6); if (T9 < .1) T9 = .1; M9 = 2.5*G1*K5; if (M9 > PIH) M9 = PIH; M9 = sin(M9); M9 = 1.+2.5*M9*sqrt(M9); if (T4 < T) goto s1680; if ((hour-T)*(T4-hour) > 0.) goto s1690; goto s1820; s1680: if ((hour-T4)*(T-hour) > 0.) goto s1820; s1690: T6 = hour+12.*(1.+sgn(T-hour))*sgn(fabs(T-hour)); 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; goto s1770; s1770: G2 = (1.+ssn/250.)*M9*sqrt(6.+58.*sqrt(G0)); G2 = G2*(1.-.1*exp((K9-24.)/3.)); G2 = G2*(1.+(1.-sgn(lat1)*sgn(lat2))*.1); G2 = G2*(1.-.1*(1.+sgn(fabs(sin(L0))-cos(L0)))); goto s1880; s1820: T6 = hour+12.*(1.+sgn(T4-hour))*sgn(fabs(T4-hour)); G8 = PI*T9/K9; U = (T4-T6)/2.; U1 = -K9/T9; G0 = C0*(G8*(exp(U1)+1.))*exp(U)/(1.+G8*G8); goto s1770; s1880: if (MUF > G2) MUF = G2; } return (MUF); } /* sgn function (like BASIC) */ double sgn(double x) { if (x > 0.) return (1.); if (x < 0.) return (-1.); return (0); }