Attribute VB_Name = "Module3"
Option Explicit

' Astrophotography utilities. This module is dependent on the
' helper functions contained in Module 1.
'   K. Fisher 5/2006 fisherka@csolutions.net

Public Function PhotogObjectArea2MPAS(ByVal objMagnitude_tmp, ByVal objType_tmp As String, ByVal objMajorAxis_arcmin_tmp, ByVal objMinorAxis_arcmin_tmp, ByVal B_mpas_tmp As Double) As String
' Accepts Magnitude and dimensions and background B_mpas.
' Returns object brightness in MPAS and contrast index.

' Valid object types are " Extended | Stellar | Stellar-CI "
'   Where object type is stellar,  objMajorAxis_arcmin_tmp and objMinorAxis_arcmin_tmp are ignored.

' Source: Covington, M. Astrophotography for Amateurs, Appendix A
'         Clark's Visual Astronomy of the Night Sky, Appendix E
'         http://clarkvision.com/visastro/appendix-e.html

'            Departs from Covington extended object model for stellar sources.
'            For a stellar source, the V-mag is the MPSAS without adjustment for area.

' Tests and usage
'    Test against objects in Clark's Appendix E catalogue
'    Clark assumes dark sky background of 24.25 MPAS
' NGC7814
' ? PhotogObjectArea2MPAS(12, 1, 5, 24.25)
' 22.3759059661737 0.749637613530523
'  Checks against Clark's results

' NGC40
' ? PhotogObjectArea2MPAS(10.5, 1, 0.7, 24.25)
' 18.7412260553693 2.20350957785228
'  Checks against Clark's results

' NGC45
' ? PhotogObjectArea2MPAS(7.8, 25, 40, 24.25)
' 23.9284809553336 0.128607617866541
'  Checks against Clark's results

' Define variables
Dim objMagnitude, objMajorAxis_arcmin, objMinorAxis_arcmin, B_mpas As Double
Dim objMajorAxis_arcsec, objMinorAxis_arcsec
Dim objB_mpas, ContrastIndex As Double
Dim strOut As String
Dim objType As String   ' Object type - Extended or Stellar
Dim ContrastIndexStrOut ' Temporary working variable for output
Dim strSpacer ' Temporary working variable for output

' Get input
objMagnitude = objMagnitude_tmp
objType = Trim(objType_tmp)
objMajorAxis_arcmin = objMajorAxis_arcmin_tmp
objMinorAxis_arcmin = objMinorAxis_arcmin_tmp
B_mpas = B_mpas_tmp

' Convert object dimensions to arcsecs
objMajorAxis_arcsec = 60# * objMajorAxis_arcmin
objMinorAxis_arcsec = 60# * objMinorAxis_arcmin

'     A) Convert objective size and v-mag to B_mpsas
'            Departs from Covington extended object model for stellar sources.
'            For a stellar source, the V-mag is the MPSAS without adjustment for area.

Select Case objType
    Case "Extended"
       ' Compute results
       objB_mpas = objMagnitude + (2.5 * (Log10((PI() / 4) * objMajorAxis_arcsec * objMinorAxis_arcsec)))
       '  Contrast index formula is from http://clarkvision.com/visastro/appendix-e.html
       ContrastIndex = -0.4 * (objB_mpas - B_mpas)
       ContrastIndexStrOut = CStr(ContrastIndex)
       strSpacer = " "
    Case "Stellar"
        objB_mpas = objMagnitude
        ContrastIndex = -0.4 * (objB_mpas - B_mpas)
        ContrastIndexStrOut = CStr(ContrastIndex)
        strSpacer = " "
    Case "Stellar-CI"
        objB_mpas = objMagnitude
        ContrastIndex = -0.4 * (objB_mpas - B_mpas)
        ContrastIndexStrOut = CStr(ContrastIndex)
        strSpacer = " "
    Case "Not applicable"      ' Should not be submitted to function, gracefully return zero
        objB_mpas = ""
        ContrastIndexStrOut = ""
        strSpacer = ""
    Case "Incomplete"      ' Should not be submitted to function, gracefully return zero
        objB_mpas = ""
        ContrastIndexStrOut = ""
        strSpacer = ""
    Case Else              ' Handle general error with graceful return of zero
        objB_mpas = ""
        ContrastIndexStrOut = ""
        strSpacer = ""
End Select

PhotogObjectArea2MPAS = Trim(CStr(objB_mpas)) & strSpacer & Trim(ContrastIndexStrOut)

End Function

Public Function PhotogStellar2MPAS(ByVal objMagnitude_tmp, ByVal B_mpas_tmp As Double) As String
' Deprecated function.  See PhotogObjectArea2MPAS, option "Stellar"
' Accepts stellar V magnitude and background B_mpas.
' Returns object brightness in MPAS and contrast index.

'  The object size is effective zero, so the returned stellar mpsas is the same as the magnitude

' Source: Covington, M. Astrophotography for Amateurs, Appendix A
'         Clark's Visual Astronomy of the Night Sky, Appendix E
'         http://clarkvision.com/visastro/appendix-e.html

' Tests and usage

' Define variables
Dim objMagnitude, objMajorAxis_arcmin, objMinorAxis_arcmin, B_mpas As Double
Dim objMajorAxis_arcsec, objMinorAxis_arcsec
Dim objB_mpas, ContrastIndex As Double
Dim strOut As String

Dim objMajorAxis_arcmin_tmp, objMinorAxis_arcmin_tmp As Double
objMajorAxis_arcmin_tmp = 0#
objMinorAxis_arcmin_tmp = 0#

' Get input
objMagnitude = objMagnitude_tmp
objMajorAxis_arcmin = objMajorAxis_arcmin_tmp
objMinorAxis_arcmin = objMinorAxis_arcmin_tmp
B_mpas = B_mpas_tmp

' Convert object dimensions to arcsecs
objMajorAxis_arcsec = 60# * objMajorAxis_arcmin
objMinorAxis_arcsec = 60# * objMinorAxis_arcmin

' Compute results
objB_mpas = objMagnitude + (2.5 * (Log10((PI() / 4) * objMajorAxis_arcsec * objMinorAxis_arcsec)))
'  Contrast index formula is from http://clarkvision.com/visastro/appendix-e.html
ContrastIndex = -0.4 * (objB_mpas - B_mpas)

PhotogStellar2MPAS = Trim(CStr(objB_mpas)) & " " & Trim(CStr(ContrastIndex))

End Function

Public Function PhotogNELMStellar2Ba(ByVal NELM_tmp As Double) As Double
' Using Schaefer relationship, converts NELM into MPSAS
'   MPAS or MPSAS = magnitudes per square arcsec squared
'   NELM is Naked-Eye Limiting Magnitude
' B_mpsas = 21.58 - 5 log(10^(1.586-NELM/5)-1) (2.0)

' Inputs: Range is limited to max of about 7.0. The limit of NELM is 6.6.

'    B_mpsas is also denoted as "Ba".  It should be distinguished from "B"
'    which is the photograph system brightness in candelas per square foot.

' Source:
' Olof Carlin, Nils. About Bradley E. Schaefer: Telescopic limiting
' Magnitudes . . . . Web page discussion of brightness in Schaefer
' (1990) and Clark (1994) at:
' http://w1.411.telia.com/~u41105032/visual/Schaefer.htm (accessed 7/2003)
' See Schaefer, Bradley E. 1990. Telescopic limiting Magnitudes. Pub. ASP 102:212-229.
' http://adsbit.harvard.edu/cgi-bin/nph-iarticle_query?bibcode=1990PASP..102..212S&db_key=AST&page_ind=0&plate_select=NO&data_type=GIF&type=SCREEN_GIF

' Test:
' ? PhotogNELMStellar2Ba(4.5)
' 18.6510697666488

' Define variables
Dim NELM As Double

' Get input
NELM = NELM_tmp

' Return Result
PhotogNELMStellar2Ba = 21.58 - (5 * Log10(10 ^ (1.586 - NELM / 5) - 1))

End Function

Public Function PhotogBa2NELMStellar(ByVal B_mpsas_tmp As Double) As Double
' Using inverse of Schaefer relationship, converts MPAS back into NELM
'   MPAS or MPSAS = magnitudes per square arcsec squared
'   NELM is Naked-Eye Limiting Magnitude
' NELM=7.93-5*log(10^(4.316-(Bmpsas/5))+1)   (1.0)
' Is inverse of Schaefer relationship.  See PhotogNELM2MPASStellar

' Outputs: Range is limited to max of about 7.0. The limit of NELM is 6.6.

' Source:
' Olof Carlin, Nils. About Bradley E. Schaefer: Telescopic limiting
' Magnitudes . . . . Web page discussion of brightness in Schaefer
' (1990) and Clark (1994) at:
' http://w1.411.telia.com/~u41105032/visual/Schaefer.htm (accessed 7/2003)

' Test:
' ? PhotogBa2NELMStellar(18.651069766488)
' 4.49999999987233

' Define variables
'    B_mpsas is also denoted as "Ba".  It should be distinguished from "B"
'    which is the photograph system brightness in candelas per square foot.
Dim Ba As Double

' Get input
Ba = B_mpsas_tmp

' Return result
PhotogBa2NELMStellar = 7.93 - (5 * Log10(10 ^ (4.316 - (Ba / 5)) + 1))

End Function

Public Function PhotogMPAS2B_candelasqft(ByVal B_mpsas_tmp As Double) As Double
' Accepts astronomical brightnesses in units of magnitudes per square arcsec
' Returns object brightness in photographic system of candelas per square foot
' B = 2.512 ^ (9.0-m")

'    B_mpsas is also denoted as "Ba".  It should be distinguished from "B"
'    which is the photograph system brightness in candelas per square foot.

' Source: Covington, M. (1998, 2d). Astrophotography for Amateurs.
'          Appendix A, at 259-261

' Tested against Covington Appendix A Tables given m" and B values
' for various planets.

' Moon (thin crescent) m" = 6.7 B = 8.0
' ?PhotogMPAS2B_candelasqft(6.7)
' 8.31850267586248

' Moon (full) m" = 3.4 B = 180
' ? PhotogMPAS2B_candelasqft(3.4)
' 173.824086762048

' Mars m" = 4.0 B = 100
' ? PhotogMPAS2B_candelasqft(4#)
'  100.022608259449

' Nebulae, faint (California) m" = 21.0 B = 0.000016 (underexpose)
' ? PhotogMPAS2B_candelasqft(21#)
' 1.58403356263317E-05

' Galaxies, bright cores m" = 18.0  B = 0.00025
' ? PhotogMPAS2B_candelasqft(18#)
' 2.51086454611934E-04

' Define variables
Dim B_mpsas As Double
' Get input
B_mpsas = B_mpsas_tmp
         
PhotogMPAS2B_candelasqft = CDbl(2.512 ^ (9# - B_mpsas))

End Function

Public Function PhotogB_candelasqft2MPAS(ByVal B_tmp As Double) As Double
' Accepts astronomical brightnesses in units of photographic system of
' candelas per square foot.  Returns object brightness in astronomical
' brightness units of magnitudes per square arcsec.

'    B_mpsas is also denoted as "Ba".  It should be distinguished from "B"
'    which is the photograph system brightness in candelas per square foot.

' Source: Covington, M. (1998, 2d). Astrophotography for Amateurs.
'          Appendix A, at 259-261
         
' Usage and test
' Moon (thin crescent) m" = 6.7 B = 8.0
' ?PhotogMPAS2B_candelasqft(6.7)
' 8.31850267586248

' Mars m" = 4.0 B = 100
' ? PhotogB_candelasqft2MPAS(100)
' 3.99878517801693

' Nebulae, faint (California) m" = 21.0 B = 0.000016 (underexpose)
' ? PhotogB_candelasqft2MPAS(0.000016)
' 20.9926131135968

' Galaxies, bright cores m" = 18.0  B = 0.00025
' ? PhotogB_candelasqft2MPAS(0.000025)
'  20.5079453201423
      
' Define variables
Dim B As Double
' Get input
B = B_tmp
         
PhotogB_candelasqft2MPAS = CDbl(9# - (1.086 * Log(B)))

End Function

Public Function PhotogExposureTimeFilmObj(ByVal objVMag_tmp, ByVal objType_tmp As String, ByVal objMajorAxis_arcmin_tmp, ByVal objMinorAxis_arcmin_tmp, ByVal Bsky_mpsas_tmp As Double, _
    ByVal dblEFR_tmp As Double, ByVal dblFilterFactor_tmp, ByVal ISOSpeed_tmp As Double, ByVal Recipority_Factor_tmp As Double) As String
' Astrophotography film exposure time for extended and stellar objects

' APPEARS TO IMPLEMENTED PROPERLY, BUT ONLY ROUGHLY AGREES WITH COVINGTON TABLES.
'   COVINGTON'S TABLES ARE EXPERIENCE ADJUSTED. NEEDS TESTING AGAINST ONLINE CALCULATORS.

' Valid object types are " Extended | Stellar | Stellar-CI "
'   Where object type is stellar,  objMajorAxis_arcmin_tmp and objMinorAxis_arcmin_tmp are ignored.

' Departs from Covington extended object model for stellar sources.  For a stellar source, the V-mag is the MPSAS without adjustment for area.

'  Receives object brightness and size information, sky info, effective focal ratio of telescope, and film recipority error
'  Returns a initial suggested exposure time in seconds and
'          an estimated sky fog time in seconds

' Source: Covington, M. Astrophotography for Amateurs, Appendix A

' Inputs
'   Object magnitude  objVMag_tmp
'   Major axis size arcmin  objMajorAxis_arcmin_tmp
'   Minor axis size arcmin objMinorAxis_arcmin_tmp
'   Sky brightness in MPSAS Bsky_mpsas_tmp
'   Effective focal ratio dblEFR_tmp
'   Film speed
'   Film's recipority factor - or Schwarzschild exponent
'        Recipority factor for typical films are
'           1.0 Digital CCD
'           0.95 Slide film
'           0.9 Newer slow films
'           0.8 older slow films
'           0.8 newer fast films
'           0.7 older fast films

' Test and usage

' Define variables
Dim objVMag, objMajorAxis_arcmin, objMinorAxis_arcmin, Bsky_mpsas As Double
Dim dblEFR, dblFilterFactor, ISOSpeed, Recipority_Factor As Double
Dim A, B, C, D, E As Double ' Working variables
Dim strOut As String
Dim M, n, P, Q, R, S As String ' Working variables
Dim Bobj_mpsas As Double          ' Object brightness in MPSAS
Dim Bobj_candelaspsqft As Double  ' Object brightness in candelas/sqft
Dim Bobj_candelas_adjusted As Double        ' Object brightness in candelas/sqft
                                            '    adjusted for the filter factor
Dim Bsky_candelaspsqft                      ' Sky brightness in candelas/sqft
Dim Bsky_candelas_adjusted                  ' Sky brightness in candelas/sqft
Dim Bsky_candelaspsqft_adjusted As Double   ' Sky brightness in candelas/sqft
                                            '    adjusted for the filter factor
Dim objType As String                       ' Object type - Extended or Stellar
                                  
                                  
' Get input
objVMag = objVMag_tmp
objType = Trim(objType_tmp)
objMajorAxis_arcmin = objMajorAxis_arcmin_tmp
objMinorAxis_arcmin = objMinorAxis_arcmin_tmp
Bsky_mpsas = Bsky_mpsas_tmp
dblEFR = dblEFR_tmp
dblFilterFactor = dblFilterFactor_tmp
ISOSpeed = ISOSpeed_tmp
Recipority_Factor = Recipority_Factor_tmp

' Find the object exposure time

' 1) Find the object B_candelaspsqft from the object V magnitude and size
'     A) Convert objective size and v-mag to B_mpsas
'            Departs from Covington extended object model for stellar sources.
'            For a stellar source, the V-mag is the MPSAS without adjustment for area.

        Bobj_mpsas = CDbl(GetCSWord(PhotogObjectArea2MPAS(objVMag, objType, objMajorAxis_arcmin, objMinorAxis_arcmin, Bsky_mpsas), 1))

' Debug.Print "Bobj_mpsas", Bobj_mpsas

'     B) Convert object B_mpsas to candelas per sq/ft
Bobj_candelaspsqft = PhotogMPAS2B_candelasqft(Bobj_mpsas)

' Debug.Print "Bobj_candelaspsqft", Bobj_candelaspsqft

' 2) Adjust B_candelaspsqft for the filter factor
Bobj_candelas_adjusted = Bobj_candelaspsqft / dblFilterFactor

' 3) Compute the initial exposure time
A = (dblEFR ^ 2) / (ISOSpeed * Bobj_candelas_adjusted)

' Debug.Print "Unadjusted object exposure time", A, A / 60#, A / 3600#
' Stop

' 4) Adjust for the initial exposure time for the reciprocity failure
B = ((A + 1#) ^ (1# / Recipority_Factor)) - 1#
' Debug.Print "Recipority adjusted exposure time", B, B / 60#, B / 3600#

' Find the sky background brightness

' 1) Find the sky brightness B_candelaspsqft from the sky Ba_mpsas
'     B) Convert object B_mpsas to candelas per sq/ft
' Debug.Print "Bsky_mpsas", Bsky_mpsas

Bsky_candelaspsqft = PhotogMPAS2B_candelasqft(Bsky_mpsas)

' Debug.Print "Bsky_candelaspsqft ", Bsky_candelaspsqft

' 2) Adjust B_candelaspsqft for the filter factor
Bsky_candelas_adjusted = Bsky_candelaspsqft / dblFilterFactor

' 3) Compute the initial exposure time
C = (dblEFR ^ 2) / (ISOSpeed * Bsky_candelas_adjusted)

' Debug.Print "Unadjusted sky exposure time", C, C / 60#, C / 3600#

' 4) Adjust for the initial exposure time for the reciprocity failure
D = ((C + 1#) ^ (1 / Recipority_Factor)) - 1#

' Debug.Print "Recipority sky adjusted exposure time", D, D / 60#, D / 3600#

'Build the output
'   Add the film exposure time to the output string
strOut = Trim(CStr(Format(B, "0.0000"))) & " " & Trim(CStr(Format(B / 60#, "0.00000"))) & " " & Trim(CStr(Format(B / 3600#, "0.000000"))) & " "
'   Add the sky fog time to the output string
strOut = strOut & Trim(CStr(Format(D, "0.0000"))) & " " & Trim(CStr(Format(D / 60#, "0.00000"))) & " " & Trim(CStr(Format(D / 3600#, "0.000000")))

' Return the result
PhotogExposureTimeFilmObj = strOut

End Function

Public Function PhotogCalibrateDigitalISOSpeed(ByVal Bsky_mpsas_tmp As Double, ByVal Exposure_time_tmp As Double, ByVal Well_Capacity_tmp As Double, ByVal Well_Exposed_tmp As Double, ByVal dblEFR_tmp As Double) As Double
' Does a down and dirty estimate of a digital camera's equivalent ISO speed
' Effective focal ratio is assumed to be 1 with the camera sitting on table and pointing to the sky.

'  Receives sky info and exposure info and returns an
'  estimte of a digital camera's ISO speed

' Source: author after HAIP at 122

' Inputs
'   Sky brightness in MPSAS Bsky_mpsas_tmp
'   Exposure time
'   Effective focal ratio
'   Well_capacity exposed and Well_capacity

' Test and usage

' Define variables
Dim Well_Capacity, Well_Exposed, Bsky_mpsas, ISO_speed, dblEFR, Well_exposure_ratio, Exposure_time As Double
Dim A, B, C, D, E As Double ' Working variables
Dim strOut As String
Dim Bsky_candelaspsqft   ' Sky brightness in candelas/sqft
                                  
' Get input
Well_Capacity = Well_Capacity_tmp
Well_Exposed = Well_Exposed_tmp
Bsky_mpsas = Bsky_mpsas_tmp
Exposure_time = Exposure_time_tmp
dblEFR = dblEFR_tmp

' Find the object exposure speed

' 1) Find the sky brigthness in photographic candelas per square foot
Bsky_candelaspsqft = PhotogMPAS2B_candelasqft(Bsky_mpsas)

' 2) Compute the well-exposure ratio
Well_exposure_ratio = Well_Exposed / Well_Capacity

' 3) Compute the adjusted time to the full exposure
' assuming a linear chip response
B = Exposure_time / Well_exposure_ratio

' 4) Estimate the ISO speed at full exposure
ISO_speed = ((dblEFR ^ 2) / (Bsky_candelaspsqft * B))

' Return result
PhotogCalibrateDigitalISOSpeed = ISO_speed

End Function


Public Function PhotogVmag2MPASStellar(ByVal Vmag_tmp As Double) As Double
' Using Schaefer relationship, converts Vmag into MPSAS
'   MPAS or MPSAS = magnitudes per square arcsec squared
'   Vmag is Johnson V filter magnitude in astronomical catalogues.

'    B_mpsas is also denoted as "Ba".  It should be distinguished from "B"
'    which is the photograph system brightness in candelas per square foot.

' Source:
' Covington, Micheal.  1999 (2d). Appendix A. Astrophotography for the amateur. p. 259.

' Test:

' Define variables
Dim Vmag As Double

' Get input
Vmag = Vmag_tmp

' Return Result
PhotogVmag2MPASStellar = 24.19 - 2.814 * Vmag + 0.3694 * Vmag ^ 2
End Function

Public Function PhotoLimitingMagQuick(ByVal F_ratio As Double) As Double
'  North's photographic limiting magnitude estimate
' Source: North, Gerald. 1997.  Advanced Amateur Astronomy. at 112

' Test data
'   None

PhotoLimitingMagQuick = 5 + (5 * Log(F_ratio))

End Function


Public Function ResolutionFilmPlusTele(ByVal Detector_lpmm_tmp As Double, ByVal Telescope_lpmm_tmp As Double) As Double
' Receives resolution of film detector in lines per mm and
' telescope in lines per mm and returns combined resolution of
' this two component system.  Is approximate because it does not include
' the MTF transfer curve.

' Source:  Covington (1999).  Astrophotography for Amateurs. at 85

' Test data - none

' Working variables
Dim Detector_lpmm, Telescope_lpmm, A As Double

' Get input
Detector_lpmm = Detector_lpmm_tmp
Telescope_lpmm = Telescope_lpmm_tmp

' Find reciporical of resolving power
A = (1# / Detector_lpmm) + (1# / Telescope_lpmm)

' Convert reciporical and return the result
ResolutionFilmPlusTele = 1# / A

End Function


Public Function PhotoRecommendFocalLength(ByVal D_mm As Double, ByVal lambda As Double, ByVal Detector_majoraxis_microns As Double) As String
' Receives an aperature and light wavelength and recommends focal length
' that optimizes the current telescope's focal length to a CCD detector size

' Berry, HAIP (2005, 2d) at 82-83.

' Dependencies:
' DiffractionDiskFWHMAngularDia_arcsec

' Test data
' Berry, HAIP (2005, 2d) at 82-83. The recommended focal length for
' a telescope of 250mm aperture viewing 600nm light, and imaging
' with a detectorat 7.5 microns, and a focal ratio of f/25 is
' 6180mm
'  Results
'  ? PhotoRecommendFocalLength(250, 600, 7.5)
'  6127.45098039215 24.5098039215686
'  6127 appears to be the correct answer, not 6180mm per HAIP p. 83

Dim A, B, C As Double
Dim D_meter, lambda_meter, Detector_majoraxis_meters As Double
Dim FWHM_arcsec, FWHM_pixel_arcsec As Double

'Convert input units to same basis - meters
 D_meter = D_mm / 1000#
 lambda_meter = (lambda * (10# ^ -9#))
 Detector_majoraxis_meters = Detector_majoraxis_microns * (10# ^ -6#)

FWHM_arcsec = DiffractionDiskFWHMAngularDia_arcsec(D_mm, lambda)

FWHM_pixel_arcsec = FWHM_arcsec / 2#

A = 206265 * Detector_majoraxis_meters / FWHM_pixel_arcsec

B = A * (10 ^ 3)   '  Recommended focal length in millimeters

C = B / D_mm '  Recommended focal ratio

PhotoRecommendFocalLength = CStr(B) & " " & CStr(C) ' Return output string

End Function


Public Function PhotoDecDrift(ByVal e_deg_tmp As Double, ByVal h_dechrs_tmp As Double, ByVal DecObj_decdeg_tmp As Double, ByVal t_min_tmp As Double) As Double
' Receives polar misalignment error and declination of target star and returns
' declination drift in arcsecs for a given period of time

' Source: Covington (1999). Astrophotography for the Amateur. At 282 per 276-279.

' Test data
'  Covington Example. 2 at 278. Worst case misalignment 1.0 deg polar alignment error,
'    object declination at celestial equator decObj = 0, time is 1 minutes,
'    or 0.25 degrees in 1 minute; object hour angle is 90 degrees (6 hours).
'    Result is 15.7 arcsecs drift in 15 minutes.
' ? PhotoDecDrift(1,6,0,1)
'  15.7071159683767

' Working variables
Dim x, y, z, A, B, C, D, f, g, Q As Double
Dim e_deg, h_dechrs, DecObj_decdeg, t_min As Double
Dim e_rads, h_rads, t_rads, ht_rads, DecObj_rads As Double

' Get input
e_deg = e_deg_tmp
h_dechrs = h_dechrs_tmp
DecObj_decdeg = DecObj_decdeg_tmp
t_min = t_min_tmp


' convert to radians
f = CDbl(e_deg)
e_rads = Deg2Rad(f)
h_rads = Deg2Rad(h_dechrs * 15#)
' Q = (h_dechrs * 15) + (t_min / 60#)
Q = (h_dechrs * 15) + ((t_min / 60#) * 15#)
ht_rads = Deg2Rad(Q)
' t_rads = Deg2Rad(t_min / 4#)
DecObj_rads = Deg2Rad(DecObj_decdeg)

' Core of computation
x = Sin(DecObj_rads) * Cos(e_rads)
y = Cos(DecObj_rads) * Sin(e_rads)
A = (Cos(ht_rads) * y) + x
B = ArcCos(A)
B = (Rad2Deg(B))  '  Convert to degrees
C = (Cos(h_rads) * y) + x
D = ArcCos(C)
D = (Rad2Deg(D))  '  Convert to degrees
z = B - D

' Discard result if too small
If (Abs(z) < 0.0001) Then
   z = 0#
Else
   ' Nothing
End If

' Return result
' Convert to arcsecs
PhotoDecDrift = z * 3600#  ' Convert to arcsecs

' Covingtons's unaltered original code required for redistribution
'  // THE COMPUTATION ITSELF
'  // Library functions and constants
'  double cos(double x)    { return Math.cos(x);  }
'  double sin(double x)    { return Math.sin(x);  }
'  double acos(double x)   { return Math.acos(x); }
'  double asin(double x)   { return Math.asin(x); }
'  double square(double x) { return x*x; }
'  double sqrt(double x)   { return Math.sqrt(x); }
'  double abs(double x)    { return Math.abs(x);  }
'  double DR = 180/Math.PI;   // degrees-to-radians factor, ~57.3
'  double decdrift(   // Declination drift, in arc-seconds, computed from:
'    double e,        //  Polar misalignment, in degrees
'    double h,        //  Hr angle of star - hr angle of false pole, degrees
'    double delta,    //  Declination of star, in degrees
'    double t)        //  Duration of exposure, in minutes of time
'  {
'    double x, y, z;                   // local variables
'    e = e/DR;                         // convert to radians
'    h = h/DR;
'    delta = delta/DR;
'    t = (t/4)/DR;
'    x = sin(delta)*cos(e);            // core of computation
'    y = cos(delta)*sin(e);
'    z = acos(cos(h+t)*y + x) - acos(cos(h)*y + x);
'    z = z * DR * 3600;                // convert to arc-seconds
'    if (abs(z) < 0.0001) z = 0.0;     // discard if too small (imprecise)
'    return z;
'  }

End Function


Public Function PhotoFieldRot(ByVal e_deg_tmp As Double, ByVal h_dechrs_tmp As Double, ByVal DecObj_decdeg_tmp As Double, ByVal t_min_tmp As Double) As Double
' Receives polar misalignment error and declination of target star and returns
' field rotation in arcsecs for a given period of time

' Source: Covington (1999). Astrophotography for the Amateur. At 282.

' Test data
'  Covington Example. 2 at 278 and 280. Worst case misalignment 1.0 deg polar alignment error,
'    object declination at celestial equator decObj = 0, time is 1 hour,
'    object hour angle is 0 degrees.  Result is 0.26 degrees of field rotation in 1 hour
' PhotoFieldRot(1, 0, 0, 60)
' 931.83684339484
' 931.83684339484 / 3600#
' 0.258843567609678

'  Covington Example. 2 at 280. Worst case misalignment 1.0 deg polar alignment error,
'    object declination at celestial equator decObj = 0, time is 1 hour or 15 degrees,
'    object hour angle is 90 degrees (6 hours).  Result is -0.034 degrees of field rotation in 1 hour
' ? PhotoFieldRot(1, 6, 0, 60)
' -122.643374273832
' ? -122.643374273831 / 3600#
' -3.40676039649531E-02

' Working variables
Dim x, y, z, w, A, B, C, D, E, f, g, H, Q As Double
Dim e_deg, h_dechrs, DecObj_decdeg, t_min As Double
Dim e_rads, h_rads, t_rads, ht_rads, DecObj_rads As Double

' Get input
e_deg = e_deg_tmp
h_dechrs = h_dechrs_tmp
DecObj_decdeg = DecObj_decdeg_tmp
t_min = t_min_tmp

' convert to radians
H = CDbl(e_deg)
e_rads = Deg2Rad(H)
h_rads = Deg2Rad(h_dechrs * 15#)
Q = (h_dechrs * 15) + ((t_min / 60#) * 15#)
ht_rads = Deg2Rad(Q)
DecObj_rads = Deg2Rad(DecObj_decdeg)

' Core of computation
'    z = asin(w*sin(h+t) / sqrt(1 - square(cos(h+t)*y + x)))
'        - asin(w*sin(h) / sqrt(1 - square(cos(h)*y + x)));

w = Sin(e_rads)
x = Sin(DecObj_rads) * Cos(e_rads)
y = Cos(DecObj_rads) * w
A = w * Sin(ht_rads)
B = 1 - ((Cos(ht_rads) * y) + x) ^ 2
C = ArcSin(A / Sqr(B))
C = Rad2Deg(C)
D = w * Sin(h_rads)
E = 1 - ((Cos(h_rads) * y) + x) ^ 2
f = ArcSin(D / Sqr(E))
f = Rad2Deg(f)
z = C - f       ' End result in degrees

' Discard result if too small
If (Abs(z) < 0.0001) Then
   z = 0#
Else
   ' Nothing
End If

' Return result
' Convert to arcsecs
PhotoFieldRot = z * 3600#  ' Convert to arcsecs

' Covingtons's unaltered original code required for redistribution
'  double rotation(   // Field rotation, in degrees
'    double e,        //  Polar misalignment, in degrees
'    double h,        //  Hr angle of star - hr angle of false pole, degrees
'    double delta,    //  Declination of star, in degrees
'    double t)        //  Duration of exposure, in minutes of time
'  {
'    double w, x, y, z;                // local variables
'    e = e/DR;                         // convert to radians
'    h = h/DR;
'    delta = delta/DR;
'    t = (t/4)/DR;
'    w = sin(e);                       // core of computation
'    x = sin(delta)*cos(e);
'    y = cos(delta)*w;
'    z = asin(w*sin(h+t) / sqrt(1 - square(cos(h+t)*y + x)))
'        - asin(w*sin(h) / sqrt(1 - square(cos(h)*y + x)));
'    z = z * DR;                       // convert to degrees
'    if (abs(z) < 0.0001) z = 0.0;     // discard if too small (imprecise)
'    return z;
'  }

End Function

Public Function PhotoFieldRotMaxPossible(ByVal e_deg_tmp As Double, ByVal t_min_tmp As Double) As Double
' Covington's maximum possible field rotation
'  Preconditions:
'    Guide star is in middle of TFOV, not off detector to one-side
'    Object declination is between +60 and -60

' Source: Covington (1999).  Astrophotography for the Amateur. At 276.

' Receives alignment error and exposure time and returns
'   an estimate of the maximum possible field rotation in arcsecs
' This function does not return a tolerlance limit of acceptable rotation

' Working variables
Dim A As Double
Dim e_deg, t_min As Double

' Get input
e_deg = e_deg_tmp
t_min = t_min_tmp

' Compute limit
A = 0.01 * e_deg * t_min

' Convert degrees to arcsecs
A = A * 3600#

' Return result
PhotoFieldRotMaxPossible = A

End Function

Public Function PhotoDecDriftMaxPossible(ByVal e_deg_tmp As Double, ByVal t_min_tmp As Double) As Double
' Covington's maximum possible acceptable declination limit
'  Preconditions:
'    Guide star is in middle of TFOV, not off detector to one-side
'    Object declination is between +80 and -80

' Receives alignment error and exposure time and returns
'   a rough estimate of maximum possible drift in arcsecs
' This function does not return a tolerance limit of acceptable drift.

' Source: Covington (1999).  Astrophotography for the Amateur. At 276.

' Working variables
Dim A As Double
Dim e_deg, t_min As Double

' Get input
e_deg = e_deg_tmp
t_min = t_min_tmp

' Compute limit
A = 15.7 * e_deg * t_min

' Return result
PhotoDecDriftMaxPossible = A

End Function



Public Function PhotoEquivalentExposure_secs(ByVal Timesecs_tmp As Double, ByVal F_ratio1_tmp As Double, ByVal F_ratio2_tmp As Double) As Double
' Receives first and new focal ratios and first exposure time
' Returns estimate of second exposure time

' Source: Sidgwick (1971 3ed).  Amateur Astronomer's Handbook. At p. 354
'  E2 = F2 ^ 2 / F1 ^ 2 * E1

' Test data
' Per Sidgwick, an f/5.6 scope exposing a plate for 60 minutes
' re-exposed at f/8 would require 122 minutes
' ? PhotoEquivalentExposure_secs(60*60,5.6,8)
'  7346.93877551021 / 60
'  122.448979591837

' Working variables
Dim Timesecs, F_ratio1, F_ratio2, A

' Get input
Timesecs = Timesecs_tmp
F_ratio1 = F_ratio1_tmp
F_ratio2 = F_ratio2_tmp

' Compute and return result
PhotoEquivalentExposure_secs = ((F_ratio2 ^ 2) / (F_ratio1 ^ 2)) * Timesecs

End Function