Find the intersection angle between two lines counter-clockwise around a circle ArcGIS 10.x

I am trying to create a tool for measuring the angle between two intersecting lines. I have created a buffer and clipped the lines to within this buffer and computed the angle between the lines using the formulas below:

dx = !X1a!-!X0a!

dy = !Y1a!-!Y0a!

dxa= !X1b!-!X0b!

dya= !Y1b!-!Y0b!

r = math.sqrt(math.pow(!dx!, 2) + math.pow(!dy!, 2))

ra = math.sqrt(math.pow(!dxa!, 2) + math.pow(!dya!, 2))

c = math.asin(abs((dxdya - dydxa))/(r*ra)) / math.pi * 180

((source How to find the angle between intersecting features in two featureclasses?, user: whuber))

This returns the angle between the two intersecting lines, but, I need the angle between the lines to always be between line1 counter-clockwise to line2. The tool will be used to calculate the angle between a road intersecting with a rail road.

Does anyone know the formula to do this?

The image below shows a sample of the data I am working with and the angle I would like to calculate. The diameter of the circle (buffer polygon) is 10 centimetres. There is a function in the math module that is much simpler to use:

deg = math.degrees(math.atan2(dy,dx))

Note the 2 on the end on tan that improves the original range from 0 - 90 degrees. This extended function will provide a 0 - 360 degree bearing for a slope (in radians). There are another pair of functions to convert between degrees and radians math.degrees() and math.radians(). No need to do that yourself.

You can then just subtract the full bearing of each line to get all the angles between two lines. The dy and dx can be negative values.

I can see that you might need to sort the end vertices to find out the direction of the lines. You do not explain why you need a particular angle. Is this because it is related to the existing line direction or do you always want the South facing angle?

All this is too fiddly to do in the expression calculator. The time has come to learn a bit of Python to put all the steps in a script and update the fields using a cursor.

#-------------------------------------------------------------------------------
# Name:        intersection_angle
# Purpose:     calculate intersection angles
#              should not be so hard!
# Author:      kimo
#
# Created:     03/08/2016
# Copyright:   (c) kimo 2016
# Licence:     Creative Commons 3.0 New Zealand
#-------------------------------------------------------------------------------
# method
#       input
#             railway line, road lines
#       output crossing points with angles, segments with bearings
# should select only those intersecting to remove redundant lines
# 1 intersect railways and roads to get intersection points
# 2 buffer intersection points, retain intersection id
# 3 intersect buffers with original rail and roads, not clip to tag id
#
# 4 for each line
#      use geometry object to find ends and hence bearing
#      write out to two dictionaries
# 5 use dictionaries to calculate angles and store back in points

import math
import arcpy

arcpy.env.workspace = "j:/temp/sample.gdb"
arcpy.env.overwriteOutput = True
# simplified set that have intersecting lines, no junctions at crossings
# single part and generally error free
inFeatures = ["road", "rail"]

# Get intersections as points and add extra fields for output
arcpy.analysis.Intersect(inFeatures, "crossing1", "", 0.01, "point")
print arcpy.GetMessages()

# but ESRI don't 'do' node topology any more
# start by deleting duplicate points
arcpy.DeleteIdentical_management(in_dataset="crossing1", fields="Shape", xy_tolerance="", z_tolerance="0")

# just select intersection points on railway lines, handle railway junctions later
arcpy.management.MakeFeatureLayer('crossing1','cross_lay')
arcpy.management.SelectLayerByLocation('cross_lay','INTERSECT','rail',1)
arcpy.management.CopyFeatures('cross_lay','crossing2')
arcpy.management.DeleteField("crossing2",['type'])

# Add a serial ID, use a cursor not field calculator if we are in python anyway
# easier to debug, trap errors, understand, just as fast
id = 1
with arcpy.da.UpdateCursor("crossing2",['id']) as cur:
for row in cur:
row = id
cur.updateRow(row)
id+=1

# Buffer crossings by 10 meters to get local angles
arcpy.analysis.Buffer("crossing2", "crossing10", "10 meters")

# Intersect original merge (not clip because ID attributes lost)
print arcpy.GetMessages()
print arcpy.GetMessages()
# each road, rail line will inherit crossing id and a type tag

# calculate bearings from ends at buffer intersections
dCross_rail = {}

with arcpy.da.UpdateCursor("roadrail_clip",['id','bearing','SHAPE@','type']) as cur:
for row in cur:
dx = row.lastPoint.X - row.firstPoint.X
dy = row.lastPoint.Y - row.firstPoint.Y
bearing = math.degrees(math.atan2(dy,dx))
# leave alone with original sense
row = bearing
cur.updateRow(row)
# to compare bearings by crossing ID
# need to collect bearings by ID and type
if row == 'road':
elif row == 'rail':
dCross_rail[row] = bearing
else:
pass
# no interest

# angles can be calculated from the two dictionaries
# if both entries exist
i = 0
dAngle = {}
for id in dCross_rail:
if dCross_rail.get(id,None) and dCross_road.get(id,None):
b_rail = dCross_rail[id] 