二地点の経緯度から方位角の求め方です。
方位角は、北が0度、東が90度、南が180度、西が270度です。

二地点の、経度(longitude)方向の距離Δlonと緯度(latitude)方向の距離Δlatがわかれば、正接関数の逆関数により方位角が計算できます。
緯度あたりの距離は、地球を真球のモデルで考えると地上のどこでも同じです。緯度一度あたりの距離は約111kmです。

一方、経度あたりの距離は、緯度によって異なります。赤道上では経度と同じですが、それ以外ではcos(lat)を乗じた値になります。極点ではゼロです。二地点間は十分近いので、cos(lat)の計算はどちらか片方の緯度を使って計算することができます。

北を正、南を負、東を正、南を負とし、Δlatを緯度による補正後とすると、方位角Azimuthは、
Azimuth = atan(Δlon/Δlat)*180/π



ただし、逆正接関数は-π/2～+π/2の範囲で値を返すので、北東側分（第一象限）しか計測できません。また、欲しい方位角は0～360度の範囲なので、適宜場合分けをして計算します。
Δlon>0, Δlat>0: Azimuth = atan(Δlon/Δlat)*180/π　第一象限
Δlon<0, Δlat>0: Azimuth = 360+atan(Δlon/Δlat)*180/π　第二象限
Δlat<0: Azimuth = 180+atan(Δlon/Δlat)*180/π　第三象限、第四象限
Cの標準ライブラリmath.hには、x,yの座標を入れて逆正接を計算するatan2関数があるので、これを使っても良いですね。Δlat≓0の場合でも、引数の無限大を考慮する必要がなくなります。

GPS受信機出力のNMEAフォーマットでは、緯度はdddmm.mmmm、経度はddmm.mmmmの形式で出力されます。1/10,000分がフォーマットの分解能になります。
180度を1/10,000分で正規化すると一億八百万なので、32bitの処理系ならばintで扱えます。math.h中の三角関数は、引数も戻り値もdoubleなので、適宜変換が必要です。
真球モデルでの緯度1/10,000分の距離は約0.1853mです。


誤差に関して考えてみます。
測位は、受信機により1秒毎におこなわれます。ある程度の分散があるので、何度も測位して平均値を求めます。
測位したA点（真の座標値である必要はありません）に対して、相対的なB点の位置を想定します。



B点でも平均値を求めますが、相対的に正しい位置(Relative B point)に対して、εの誤差を持つと考えられます。図では、経度方向に10、緯度方向に3の誤差がある例です。

AB間の距離が100、Aに対するBの方位角が45度（北東）であるとします。真のΔlonは100/SIN(π/4)、Δlatは100/COS(π/4)、共に100/SQRT(2)です。
上記の誤差がある場合、
Δlon = 100/SQRT(2) + 10
Δlat = 100/SQRT(2) + 3
として測位されます。このときのAzimuthは、前項の式に当てはめると、計算で得られたAzimuthは約47.60度となり、2.6度の誤差になります。方位角0～45では経度方向、45～90では緯度方向の座標誤差がより支配的になります。
いくつかのAzimuthの値、AB間距離に対する座標誤差（%）で、Azimuth誤差（分）をシミュレーションした表が以下です。



太陽の視直径（視野角）は約30分です。望遠鏡や望遠カメラに70%の視野で捉える前提であれば、30分に対して(1-0.7)/2以内の誤差、即ち4.5分（0.075度）以内である必要があります。
シミュレーション結果から、A、B点の測位相対誤差は、おおよそ0.1%以下でなければならないことがわかります。
一方、受信機のNMEAフォーマットの分解能が約0.1853mなので、AB間距離は少なくとも200m必要になってしまいます。近いAB間距離で測位をおこなうには、別のフォーマットで出力される高精度の受信機が必要ということになります。
あまり現実的な条件ではなさそうです。日食撮影中に、カメラ方向の微調整をする方向で、精度を緩和した方が良さそうです。