eddy/BTA_utils
1
0
Fork 0
mirror of https://github.com/eddyem/BTA_utils.git synced 2025-12-06 02:35:13 +03:00
Code Issues Packages Projects Releases Wiki Activity
BTA_utils/bta_print_header/bta_print.c

481 lines
16 KiB
C
Raw Blame History

/*
* This file is part of the btaprinthdr project.
* Copyright 2022 Edward V. Emelianov <edward.emelianoff@gmail.com>, Vladimir S. Shergin <vsher@sao.ru>.
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <time.h>
#include <usefull_macros.h>
#include <erfa.h>
#include <erfam.h>
#include "bta_print.h"
#include "bta_shdata.h"
#include "bta_site.h"
// rad to time sec
#ifndef ERFA_DR2S
#define ERFA_DR2S 1.3750987083139757010431557155385240879777313391975e4
#endif
// seconds of time in full circle
#define S24 (24.*3600.)
const double longitude = TELLONG * 3600.;
const double Fi = TELLAT * 3600.;
const double cos_fi=0.7235272793;
const double sin_fi=0.6902957888;
/*
const double cos_fi = 0.723527277857; // Cos of SAO latitude
const double sin_fi = 0.690295790366; // Sin --- "" -----
*/
/**
* convert apparent coordinates (nowadays) to mean (JD2000)
* appRA, appDecl in seconds
* r, d in seconds
*/
void calc_mean(double appRA, double appDecl, double *r, double *dc){
double ra=0., dec=0., utc1, utc2, tai1, tai2, tt1, tt2, fd, eo, ri;
int y, m, d, H, M;
DBG("appRa: %g'', appDecl'': %g", appRA, appDecl);
appRA *= ERFA_DS2R;
appDecl *= ERFA_DAS2R;
#define ERFA(f, ...) do{if(f(__VA_ARGS__)){WARNX("Error in " #f); goto rtn;}}while(0)
// 1. convert system JDate to UTC
ERFA(eraJd2cal, JDate, 0., &y, &m, &d, &fd);
fd *= 24.;
H = (int)fd;
fd = (fd - H)*60.;
M = (int)fd;
fd = (fd - M)*60.;
ERFA(eraDtf2d, "UTC", y, m, d, H, M, fd, &utc1, &utc2);
ERFA(eraUtctai, utc1, utc2, &tai1, &tai2);
ERFA(eraTaitt, tai1, tai2, &tt1, &tt2);
eraAtic13(appRA, appDecl, tt1, tt2, &ri, &dec, &eo);
ra = eraAnp(ri + eo);
ra *= ERFA_DR2S;
dec *= ERFA_DR2AS;
#undef ERFA
rtn:
if(r) *r = ra;
if(dc) *dc = dec;
}
static char buf[1024];
char *time_asc(double t){
int h, m;
double s;
h = (int)(t/3600.);
m = (int)((t - (double)h*3600.)/60.);
s = t - (double)h*3600. - (double)m*60.;
h %= 24;
if(s>59) s=59;
sprintf(buf, "%d:%02d:%04.1f", h,m,s);
return buf;
}
char *angle_asc(double a){
char s;
int d, min;
double sec;
if (a >= 0.)
s = '+';
else{
s = '-'; a = -a;
}
d = (int)(a/3600.);
min = (int)((a - (double)d*3600.)/60.);
sec = a - (double)d*3600. - (double)min*60.;
d %= 360;
if(sec>59.9) sec=59.9;
sprintf(buf, "%c%d:%02d:%04.1f", s,d,min,sec);
return buf;
}
/**
* @brief printhdr - write FITS record into output file
* @param fd - fd to write
* @param key - key
* @param val - value
* @param cmnt - comment
* @return 0 if all OK
*/
static int printhdr(int fd, const char *key, const char *val, const char *cmnt){
char tmp[81];
char tk[9];
if(strlen(key) > 8){
strncpy(tk, key, 8);
key = tk;
}
if(cmnt){
snprintf(tmp, 80, "%-8s= %-21s / %s", key, val, cmnt);
}else{
snprintf(tmp, 80, "%-8s= %s", key, val);
}
size_t l = strlen(tmp);
tmp[l] = '\n';
++l;
if(write(fd, tmp, l) != (ssize_t)l){
WARN("write()");
return 1;
}
return 0;
}
static void calc_AZ(double alpha, double delta, double stime, double *az, double *zd){
double ha = (stime - alpha) * 15.;
if(ha < 0.) ha += ERFA_TURNAS;
ha *= ERFA_DAS2R;
double A, alt;
eraHd2ae(ha, delta*ERFA_DAS2R, TELLAT*ERFA_DD2R, &A, &alt);
if(az) *az = (A - M_PI) * ERFA_DR2AS; // BTA counts Az from south, not north
if(zd) *zd = (M_PI_2 - alt) * ERFA_DR2AS;
}
extern void calc_airmass(
// in
double daynum, // Day number from beginning of year, 0 = midnight Jan 1
double relhumid,// Relative humidity in percent
double p0, // local pressure in mmHg
double t0, // temperature in kelvins
double z, // zenith distance in degr
// out
double *am, // AIRMASS
double *acd, // column density
double *wcd, // water vapor column density
double *wam // water vapor airmass
);
#if 0
//double coeff[8] = { -97.4, 13.0, -11.7, -3.6, -6.2, -279.9, -33.3, 8.2};
void pos_corr (double Delta, double A, double Z, double P,
double *pcA, double *pcZ) {
double ar, dr, zr, pr;
dr = Delta;
ar = A;
zr = Z;
pr = P;
double *coeff = (double*)sdt->pc_coeff;
*pcA = coeff[0] + coeff[1]/tan(zr) + coeff[2]/sin(zr) -
coeff[3]*sin(ar)/tan(zr) + coeff[4]*cos(dr)*cos(pr)/sin(zr);
*pcZ = coeff[5] + coeff[6]*sin(zr) + coeff[7]*cos(zr) +
coeff[3]*cos(ar) + coeff[4]*cos_fi*sin(ar);
}
void calc_AZP(double alpha, double delta, double stime,
double *az, double *zd, double *pa){
double sin_t,cos_t, sin_d,cos_d, sin_z,cos_z, sin_p,cos_p;
double t, d, z, a, p , x, y;
t = (stime - alpha) * 15.;
red("st=%g, alpha=%g, delta=%g, t=%g\n", stime/3600., alpha/3600., delta/3600., t/3600./15.);
if (t < 0.)
t += ERFA_TURNAS;
t *= ERFA_DAS2R;
d = delta * ERFA_DAS2R;
sincos(t, &sin_t, &cos_t);
sincos(d, &sin_d, &cos_d);
cos_z = cos_fi * cos_d * cos_t + sin_fi * sin_d;
z = acos(cos_z);
sin_z = sqrt(1. - cos_z*cos_z);
y = cos_d * sin_t;
x = cos_d * sin_fi * cos_t - cos_fi * sin_d;
a = atan2(y, x);
sin_p = sin_t * cos_fi / sin_z;
cos_p = (sin_fi * cos_d - sin_d * cos_fi * cos_t) / sin_z;
p = atan2(sin_p, cos_p);
if (p < 0.0)
p += 2.0*M_PI;
double pca, pcz;
pos_corr(d, a, z, p, &pca, &pcz);
red("pca=%g, pcz=%g; tel_cor_A=%g, tel_cor_Z=%g\n", pca, pcz, tel_cor_A, tel_cor_Z);
*zd = z * ERFA_DR2AS;
*az = a * ERFA_DR2AS;
*pa = p * ERFA_DR2AS;
}
#endif
#define WRHDR(k, v, c) do{if(printhdr(fd, k, v, c)){goto returning;}}while(0)
int print_header(_U_ const char *path){
int ret = FALSE;
char val[23], comment[71];
#define COMMENT(...) do{snprintf(comment, 70, __VA_ARGS__);}while(0)
#define VAL(fmt, x) do{snprintf(val, 22, fmt, x);}while(0)
#define VALD(x) VAL("%.10f", x)
#define VALS(x) VAL("'%s'", x)
int l = strlen(path) + 7;
char *aname = MALLOC(char, l);
snprintf(aname, l, "%sXXXXXX", path);
int fd = mkstemp(aname);
if(fd < 0){
WARN("Can't write header file, mkstemp()");
FREE(aname);
return FALSE;
}
fchmod(fd, 0644);
WRHDR("TELESCOP", "'BTA 6m telescope'", "Telescope name");
WRHDR("ORIGIN", "'SAO RAS, Russia'", "Organization responsible for the data");
VALD(TELLAT);
COMMENT("Telescope lattitude (degr): %s", angle_asc(TELLAT*3600.));
WRHDR("SITELAT", val, comment);
VALD(TELLONG);
COMMENT("Telescope longitude (degr): %s", angle_asc(TELLONG*3600.));
WRHDR("SITELONG", val, comment);
VAL("%.1f", TELALT);
WRHDR("SITEALT", val, "Telescope altitude (m)");
double sidtm = S_time-EE_time;
if(sidtm < 0.) sidtm += S24;
else if(sidtm > S24){
int x = (int)(sidtm / S24);
sidtm -= S24 * (double)x;
}
COMMENT("Sidereal time, seconds: %s", time_asc(sidtm));
VALD(sidtm);
WRHDR("ST", val, comment);
COMMENT("Universal time, seconds: %s", time_asc(M_time));
VALD(M_time);
WRHDR("UT", val, comment);
VALD(JDate);
WRHDR("JD", val, "Julian date");
{time_t t_now = time(NULL);
struct tm *tm_ut = gmtime(&t_now);
double l = ERFA_DTY; int y = 1900 + tm_ut->tm_year;
if((0 == y%400) || ((0 == y%4)&&(0 != y%100))) l += 1.;
double dtmp = y + (double)tm_ut->tm_yday / l;
VALD(dtmp);
WRHDR("EQUINOX", val, "Coordinates epoch");}
switch(Tel_Focus){
default:
case Prime:
VALS("Prime"); break;
case Nasmyth1:
VALS("Nasmyth1"); break;
case Nasmyth2:
VALS("Nasmyth2"); break;
}
WRHDR("FOCUS", val, "Observation focus");
char *telmode = "Undefined";
if(Tel_Hardware == Hard_Off) telmode = "Off";
else if(Tel_Mode != Automatic) telmode = "Manual";
else{
switch (Sys_Mode){
default:
case SysStop : telmode = "Stopping"; break;
case SysWait : telmode = "Waiting"; break;
case SysPointAZ :
case SysPointAD : telmode = "Pointing"; break;
case SysTrkStop :
case SysTrkStart:
case SysTrkMove :
case SysTrkSeek : telmode = "Seeking"; break;
case SysTrkOk : telmode = "Tracking"; break;
case SysTrkCorr : telmode = "Correction";break;
case SysTest : telmode = "Testing"; break;
}
}
VALS(telmode);
WRHDR("TELMODE", val, "Telescope working mode");
VAL("%.2f", val_F);
WRHDR("VAL_F", val, "Focus value of telescope (mm)");
double a2000, d2000;
calc_mean(InpAlpha, InpDelta, &a2000, &d2000);
VALD(a2000 * 15. / 3600.);
COMMENT("Input R.A. for J2000 (deg): %s", time_asc(a2000));
WRHDR("RA_INP0", val, comment);
VALD(d2000 / 3600.);
COMMENT("Input Decl. for J2000 (deg): %s", angle_asc(d2000));
WRHDR("DEC_INP0", val, comment);
calc_mean(CurAlpha, CurDelta, &a2000, &d2000);
VALD(a2000 * 15. / 3600.);
COMMENT("Telescope R.A. for J2000 (deg): %s", time_asc(a2000));
WRHDR("RA_0", val, comment);
VALD(d2000 / 3600.);
COMMENT("Telescope Decl. for J2000 (deg): %s", angle_asc(d2000));
WRHDR("DEC_0", val, comment);
#define RA(ra, dec, text, pref) do{VALD(ra*15./3600.); COMMENT(text " R.A. (degr): %s", time_asc(ra)); WRHDR("RA" pref, val, comment); \
VALD(dec/3600.); COMMENT(text " Decl (degr): %s", angle_asc(dec)); WRHDR("DEC" pref, val, comment);}while(0)
RA(InpAlpha, InpDelta, "Input", "_INP");
RA(CurAlpha, CurDelta, "Current object", "_OBJ");
RA(SrcAlpha, SrcDelta, "Source", "_SRC");
RA(val_Alp, val_Del, "Telescope", "");
#undef RA
#define AZ(a, z, text, pref) do{VALD(a/3600.); COMMENT(text " Az (degr): %s", angle_asc(a)); WRHDR("A" pref, val, comment); \
VALD(z/3600.); COMMENT(text " ZD (degr): %s", angle_asc(z)); WRHDR("Z" pref, val, comment);}while(0)
AZ(InpAzim, InpZdist, "Input", "_INP");
AZ(tag_A, tag_Z, "Current object", "_OBJ");
AZ(val_A, val_Z, "Telescope", "");
#undef AZ
/*
VALD(tag_A/3600.);
COMMENT("Target Az (degr): %s", angle_asc(tag_A));
WRHDR("TARG_A", val, comment);
VALD(tag_Z/3600.);
COMMENT("Target ZD (degr): %s", angle_asc(tag_Z));
WRHDR("TARG_Z", val, comment);
VALD(val_A/3600.);
COMMENT("Telescope Az (degr): %s", angle_asc(val_A));
WRHDR("A", val, comment);
VALD(val_Z/3600.);
COMMENT("Telescope ZD (degr): %s", angle_asc(val_Z));
WRHDR("Z", val, comment);
*/
double t = sidtm * ERFA_DS2R - val_Alp * ERFA_DS2R;
if(t < 0) t += 2*M_PI;
double P = ERFA_DR2D * eraHd2pa(t, val_Del * ERFA_DAS2R, TELLAT * ERFA_DD2R);
if(P < 0.) P += 360.;
VALD(P);
COMMENT("Parallactic angle (degr): %s", angle_asc(P*3600.));
WRHDR("PARANGLE", val, comment);
VALD(tag_P/3600.);
COMMENT("Target par. angle (degr): %s", angle_asc(tag_P));
WRHDR("TAGANGLE", val, comment);
VALD(val_P/3600.);
COMMENT("Current P2 rot. angle (degr): %s", angle_asc(val_P));
WRHDR("ROTANGLE", val, comment);
VALD(val_D/3600.);
COMMENT("Dome Az (degr): %s", angle_asc(val_D));
WRHDR("DOME_A", val, comment);
if(Sys_Mode==SysTrkSeek||Sys_Mode==SysTrkOk||Sys_Mode==SysTrkCorr||Sys_Mode==SysTrkStart||Sys_Mode==SysTrkMove){
double curA,curZ,srcA,srcZ;
double corAlp,corDel,corA,corZ;
corAlp = CurAlpha-SrcAlpha;
corDel = CurDelta-SrcDelta;
if(corAlp > 23*3600.) corAlp -= 24*3600.;
if(corAlp < -23*3600.) corAlp += 24*3600.;
calc_AZ(SrcAlpha, SrcDelta, S_time, &srcA, &srcZ);
calc_AZ(CurAlpha, CurDelta, S_time, &curA, &curZ);
corA = curA-srcA;
corZ = curZ-srcZ;
VALD(corAlp);
WRHDR("RACORR", val, "RA correction (current - source)");
VALD(corDel);
WRHDR("DECCORR", val, "DEC correction (current - source)");
VALD(corA);
WRHDR("ACORR", val, "A correction (current - source)");
VALD(corZ);
WRHDR("ZCORR", val, "Z correction (current - source)");
}
/*
double az,zd;//,pa;
calc_AZ(val_Alp, val_Del, sidtm, &az, &zd);
green("calc_AZ: AZ=%g, ZD=%g\n", az/3600., zd/3600.);
calc_AZx(val_Alp, val_Del, sidtm, &az, &zd);
green("calc_AZx: AZ=%g, ZD=%g\n", az/3600., zd/3600.);
calc_AZP(val_Alp, val_Del, sidtm, &az, &zd, &pa);
green("AZ=%g, ZD=%g, PA=%g; sidtm=%g\n", az/3600.,zd/3600.,pa/3600., sidtm/3600.);
*/
if(UsePCorr) VALS("true"); else VALS("false");
WRHDR("USEPCORR", val, "P.corr.sys.: K0..K7 (real = measured - PCS)");
if(UsePCorr){
char kx[3];
const char *descr[8] = {
"A0 - PCS Azimuth zero",
"L - PCS horizontal axe inclination",
"k - PCS collimation error",
"F - PCS lattitude error of vert. axe",
"dS - PCS time error",
"Z0 - PCS zenith zero",
"d - PCS tube bend [sin(Z)]",
"d1 - PCS tube bend [cos(Z)]"
};
for(int i = 0; i < 8; ++i){
VAL("%.2f", PosCor_Coeff[i]);
snprintf(kx, 3, "K%d", i);
WRHDR(kx, val, descr[i]);
}
#if 0
// SKN:
// dA = K0 + K1/tg(Z) + K2/sin(Z) - K3*sin(A)/tg(Z) + K4 *cos(delta)*cos(P)/sin(Z)
// dZ = K5 + K6*sin(Z) + K7*cos(Z) + K3*cos(A) + K4*cos(phi)*sin(A)
double A = val_A * ERFA_DAS2R, Z = val_Z * ERFA_DAS2R;
double sinZ, cosZ, sinA, cosA, tgZ = tan(Z);
sincos(Z, &sinZ, &cosZ);
sincos(A, &sinA, &cosA);
double dA = PosCor_Coeff[0] + PosCor_Coeff[1]/tgZ + PosCor_Coeff[2]/sinZ -
PosCor_Coeff[3]*sinA/tgZ +
PosCor_Coeff[4]*cos(CurDelta*ERFA_DAS2R)*cos(P*ERFA_DD2R)/sinZ;
double dZ = PosCor_Coeff[5] + PosCor_Coeff[6]*sinZ + PosCor_Coeff[7]*cosZ +
PosCor_Coeff[3]*cosA + PosCor_Coeff[4]*cos_fi*sinA;
red("dA=%g, dZ=%g; tel_cor_A=%g, tel_cor_Z=%g\n", dA, dZ, tel_cor_A, tel_cor_Z);
#endif
VAL("%.1f", pos_cor_A);
WRHDR("PCSDA_O", val, "Pointing correction for object A (arcsec)");
VAL("%.1f", pos_cor_Z);
WRHDR("PCSDZ_O", val, "Pointing correction for object Z (arcsec)");
VAL("%.1f", tel_cor_A);
WRHDR("PCSDA_T", val, "Pointing correction for telescope A (arcsec)");
VAL("%.1f", tel_cor_Z);
WRHDR("PCSDZ_T", val, "Pointing correction for telescope Z (arcsec)");
}
VAL("%.1f", refract_Z);
WRHDR("REFR_O", val, "Refraction for object position (arcsec)");
VAL("%.1f", tel_ref_Z);
WRHDR("REFR_T", val, "Refraction for telescope position (arcsec)");
{double refa, refb, tz = tan(val_Z*ERFA_DAS2R);
eraRefco(val_B*1.33322, val_T1, val_Hmd/100., 0.45, &refa, &refb);
double dz = (refa*tz + refb*tz*tz*tz)*ERFA_DR2AS;
VAL("%.2f", dz);
WRHDR("REFR_T_E", val, "RREFR_T by ERFA eraRefco()");}
#define T(hdr, t, text) do{VAL("%.1f", t); COMMENT(text " temperature (degC)"); WRHDR(hdr, val, comment);}while(0)
T("OUTTEMP", val_T1, "Outern");
T("DOMETEMP", val_T2, "In-dome");
T("MIRRTEMP", val_T3, "Mirror");
#undef T
VAL("%.1f", val_B);
WRHDR("PRESSURE", val, "Atm. pressure (mmHg)");
VAL("%.1f", val_Wnd);
WRHDR("WIND", val, "Wind speed (m/s)");
VAL("%.1f", val_Hmd);
WRHDR("HUMIDITY", val, "Relative humidity (%)");
/*
* Airmass calculation
* by Reed D. Meyer
*/
{double am, acd, wcd, wam;
time_t t_now = time(NULL);
struct tm *tm_loc;
tm_loc = localtime(&t_now);
calc_airmass(tm_loc->tm_yday, val_Hmd, val_B, val_T1+273.15, val_Z/3600., &am, &acd, &wcd, &wam);
VALD(am);
WRHDR("AIRMASS", val, "Air mass by Reed D. Meyer");
VALD(wam);
WRHDR("WVAM", val, "Water vapour air mass by Reed D. Meyer");
VALD(acd);
WRHDR("ATMDENS", val, "Atm. column density by Reed D. Meyer (g/cm^2)");
VALD(wcd);
WRHDR("WVDENS", val, "WV column density by Reed D. Meyer (g/cm^2)");}
ret = TRUE;
returning:
close(fd);
rename(aname, path);
FREE(aname);
return ret;
}
Reference in New Issue View Git Blame Copy Permalink