The GBM model can be written as:

$$ \delta S_t= \mu S_t \delta t+\sigma S_t\delta t $$

The above is short-hand for the following SDE:

$$ S(t)=S(0)+\int^{t}_{0}\mu S(h)dh+\int^{t}_{0}\sigma S(h)dW(h) $$ 

Solving the above SDE yields an expression that you implemented in your code:

$$ S(t)=S_0exp\left((\mu-0.5 \sigma^2)t+\sigma \sqrt{t}  Z\right) $$

The Black-Scholes formula can be derived directly by applying the option pay-off to the above solution of the SDE (below I use the real-world measure for simplicity*, see asterix note further below in the text for more details):

$$ Call(t_0)=e^{-rt}\mathbb{E}\left[ (S_t-K)I_{ \left( S_t>K \right) } \right] = \\ = e^{-rt}\mathbb{E}\left[ S_tI_{ \left( S_t>K \right) }-KI_{ \left( S_t>K \right) } \right]=\\=e^{-rt}\mathbb{E}\left[ S_tI_{ \left( S_t>K \right) }\right]-e^{-rt}K\mathbb{E}\left[ I_{ \left( S_t>K \right) }\right]$$

Focusing on the second term:

$$ e^{-rt}K\mathbb{E}\left[ I_{ \left( S_t>K \right) }\right] = e^{-rt}K\mathbb{P}\left( S_t>K \right) = \\ = e^{-rt}K\mathbb{P}\left( S_0 exp\left((\mu-0.5 \sigma^2)t+\sigma \sqrt{t}  Z\right)>K \right) = \\ = e^{-rt}K\mathbb{P}\left( (\mu-0.5 \sigma^2)t+\sigma \sqrt{t}  Z>ln \left(\frac{K}{S_0} \right) \right) = \\ = e^{-rt}K\mathbb{P}\left( Z>\frac{ln \left(\frac{K}{S_0} \right) -\mu t + 0.5\sigma^2 t}{\sigma \sqrt{t} } \right) = \\ = e^{-rt}K\mathbb{P}\left( Z> (-1)\frac{ln \left(\frac{S_0}{K} \right) +\mu t - 0.5\sigma^2 t}{\sigma \sqrt{t} } \right) = \\ = e^{-rt}K\mathbb{P}\left( Z \leq \frac{ln \left(\frac{S_0}{K} \right) +\mu t - 0.5\sigma^2 t}{\sigma \sqrt{t} } \right) = \\ = e^{-rt}KN(d_2)  $$

The first term $e^{-rt}\mathbb{E}\left[ S_tI_{ \left( S_t>K \right) }\right]$ requires a tiny little bit more work to evaluate, but using a similar technique this term comes out as $S_0N(d_1)$.

**So what this tedious usage of formulas was meant to demonstrate is that the Black-Scholes formula can be shown to be a direct consequence of the GBM model for the underlying stock price: therefore this answers your first and second questions:** 

(i) Yes, the mu and sigma in both models are identical, because the BS formula is based on the GBM model 

(ii) Yes, both models need to be consistent with one another in terms of units of time.

***Word of warning**: there is one additional step that needs to be performed when using the GBM model for pricing options: you should switch from the real world probability measure to the risk-neutral measure. In practical terms, it means that your drift $\mu$ needs to be replaced with drift $r$, where $r$ should be the "risk-free" rate corresponding to the option maturity. If you don't have access to the entire OIS curve for USD, then I would just take the FED funds rate as a proxy for $r$ (right now, the FED funds rate is 0.25%).

You should also use implied volatility to price the option. But using historical volatility (as you do in your code) as a proxy is ok if you just want to experiment.

Your **third question**: if you want to price the option by Monte-Carlo (i.e. simulating stock price first, then taking expectation of the option pay-off at maturity), you need to run "n" simulations (i.e. loops). But because you know the analytical solution to the GBM model as shown above and you can plug this directly into the option pay-off and **analytically** compute the option price that way, you don't actually need to run an MC simulation. You can just price the option directly via the B-S formula. 

It's basically up to you if you want to evaluate the *expectation* in the Option pay-off formula via Monte-Carlo or analytically (which leads to the BS formula directly). Obviously, analytical evaluation is more accurate than numerical approximation. Running a numerical simulation on a problem which you know how to solve analytically is a bit like hiding your own Easter eggs and then searching for them.

PS: last but not least, you should not use 365 days, but rather 260 days per year (because there are only roughly 260 trading days in a calendar year).