************************************************************;
* Program Name: Meteor Radiant Determination from an Image *;
* Last Updated: 09/11/19 by David Oesper                   *;
************************************************************;
data stars;
   retain star1 ra1 dc1 star2 ra2 dc2 b1 b2;
   input #2 star1 $  1-14
	    rh1     16-17
	    rm1     19-20
	    rs1     22-26
            dd1     28-30
            dm1     32-33
            ds1     35-38
            ;
   input #3 star2 $  1-14
	    rh2     16-17
	    rm2     19-20
	    rs2     22-26
            dd2     28-30
            dm2     32-33
            ds2     35-38
            ;
   ra1 = (rh1 + rm1/60 + rs1/3600)*(constant('pi')/12);
   dc1 = sign(dd1)*(abs(dd1) + dm1/60 + ds1/3600)*(constant('pi')/180);
   ra2 = (rh2 + rm2/60 + rs2/3600)*(constant('pi')/12);
   dc2 = sign(dd2)*(abs(dd2) + dm2/60 + ds2/3600)*(constant('pi')/180);
   b1 = (180/constant('pi'))*atan((cos(dc2)*sin(ra2-ra1))/(cos(dc1)*sin(dc2)-sin(dc1)*cos(dc2)*cos(ra2-ra1)));
   if b1 < 0 then b1 = b1 + 360;
   drop rh1 rm1 rs1 dd1 dm1 ds1 rh2 rm2 rs2 dd2 dm2 ds2;
datalines;
....,....1....,....2....,....3....,....4
3UC 148-239423 17 56 38.42 -16 23 27.1
3UC 147-243087 17 56 31.96 -16 30 40.0
run;

proc print data=stars;
   title 'stars';
run;

data radiants;
   input @1  rad    $3.
         @5  dt1 date9.
         @15 r1      3.
         @19 d1      3.
         @23 dt2 date9.
         @33 r2      3.
         @37 d2      3.
         @41 dtm date9.
         ;
   rar = (constant('pi')/180)*(r1 + (r2-r1)*((dtm-dt1)/(dt2-dt1)));
   dcr = (constant('pi')/180)*(d1 + (d2-d1)*((dtm-dt1)/(dt2-dt1)));
datalines;
ANT 15aug2019 335 -08 20aug2019 340 -07 16aug2019
PER 15aug2019 051 +58 20aug2019 057 +58 16aug2019
SDA 15aug2019 352 -12 20aug2019 356 -11 16aug2019
KCG 15aug2019 285 +59 20aug2019 286 +59 16aug2019
run;

proc print data=radiants;
   title 'radiants';
run;

proc sql;
   create table combine as
      select r.*
            ,s.*
      from radiants as r
          ,stars    as s
;quit;
run;

proc print data=combine;
   title 'combine';
run;

data results;
   set combine;
   br1 = (180/constant('pi'))*atan((cos(mean(dc1,dc2))*sin(mean(ra1,ra2)-rar))/(cos(dcr)*sin(mean(dc1,dc2))-sin(dcr)*cos(mean(dc1,dc2))*cos(mean(ra1,ra2)-rar)));
   if br1 < 0 then br1 = br1 + 360;
   df = abs(b1-br1);
run;

proc sort data=results;
   by df;
run;

proc print data=results;
   title 'results';
run;

proc print data=results label noobs;
   label rad = 'Meteor Shower'
         b1  = 'Field Bearing'
         br1 = 'Radiant Bearing'
         df  = 'Difference'; 
   title 'Meteor Shower Radiant Match';
   var rad b1 br1 df;
run;