import numpy as np
import scipy.stats as stats

def call(S, K, T, r, div):
    d1 = (np.log(S / K) + (r - div + (vol ** 2) / 2) * T) \
         / (vol * np.sqrt(T))
    d2 = (np.log(S / K) + (r - div - (vol ** 2) / 2) * T) \
         / (vol * np.sqrt(T))
    price = (S * np.exp(-div * T) * stats.norm.cdf(d1, 0.0, 1.0)
             - K * np.exp(-r * T) * stats.norm.cdf(d2, 0.0, 1.0))
    return price

def put(S, K, T, r, div):
    d1 = (np.log(S / K) + (r - div + (vol ** 2) / 2) * T) \
         / (vol * np.sqrt(T))
    d2 = (np.log(S / K) + (r - div - (vol ** 2) / 2) * T) \
         / (vol * np.sqrt(T))
    price = (K * np.exp(-r * T) * stats.norm.cdf(-d2, 0.00, 1.0)
             - S * np.exp(-div * T) * stats.norm.cdf(-d1, 0.0, 1.0))
    return price

option_type = input("Price a C or a P?")
market_price = float(input("Input market price: "))
S = float(input("Input spot: "))
K = float(input("Input strike: "))
T = float(input("Input time to maturity: "))
r = float(input("Interest rate: "))
div = float(input("Input dividend yield: "))
vol = 0
implied_price = 0

# Add more decimals to increase accuracy
while (market_price - implied_price) > 0.0001:
    vol += 0.0001       # Add more decimals to increase accuracy
    if option_type == 'C':
        implied_price = call(S, K, T, r, div)
    elif option_type == 'P':
        implied_price = put(S, K, T, r, div)
else:
    print("Implied vol is :" + str(vol))