(* PROTOTYP kodu, normalnie zostaje wyrzucony z wersji końcowej. Właściwe funkcje dalej. *)

H=10

V0x=1

V0y=1

V0z=7

g=10

m=1

k=0.6

V0=Sqrt[V0x^2+V0y^2+V0z^2]

rozw=NDSolve[{m x''[t]==-k x'[t],m y''[t]==-k y'[t],m z''[t]==-m g -k z''[t],x[0]==0,y[0]==0,z[0]==H,x'[0]==V0x,y'[0]==V0y,z'[0]==V0z},{x,y,z},{t,0,2*V0^2/g*Sqrt[1+2 H g/V0^2]}]//First

Z=z/.rozw

Plot[Z[t],{t,0,10}]

tZ=FindRoot[Z,{V0^2/g*Sqrt[1+2 H g/V0^2]}]//First

X=x/.rozw;
Y=y/.rozw;

{X[tZ],Y[tZ]}

(* KONIEC PROTOTYPU *)

(* 
Funkcja obliczająca współrzędne punktu, na który spada pocisk wystrzelony z prędkością o składowych V0x, V0y, V0z, z wysokości H, w polu grawitatyjnym o natężeniu g. Masa pocisku wynosi m, a współcznnik oporu powietrza k.
*)
wspolrzedne[V0x_,V0y_,V0z_,H_,g_,m_,k_]:=Module[{V0,czasLotu,rozwiazanierr,X,Y,Z,tZ},
V0=Sqrt[V0x^2+V0y^2+V0z^2];
czasLotu=V0^2/g*Sqrt[1+2 H g/V0^2];(* szacowany czas lotu pocisku *)
(* Polecenie niżej to numeryczne rozwiązanie równań Newtona *)
rozwiazanierr=NDSolve[{
m x''[t]==-k x'[t],m y''[t]==-k y'[t],m z''[t]==-m g -k z''[t],(* równania ruchu m a == F, a=x''[t] *)
x[0]==0,y[0]==0,z[0]==H,(* położenie w chwili t=0 *)
x'[0]==V0x,y'[0]==V0y,z'[0]==V0z},(* prędkości w chwili t=0 *)
{x,y,z},(* trzy szukane funkcje *)
{t,0,5*czasLotu} (* maksymalny czas ,,lotu'' *)
]//First; (* First jest potrzebne, bo wynik ma postać listy (jednoelementowej) *)
Z=z/.rozwiazanierr;(* funkcja opisująca zależność wysokości pocisku od czasu *)
tZ=FindRoot[Z,{czasLotu}]//First; (* szukamy numerycznie momentu ,,upadku'', czyli rozwiązujemy równanie Z[t]==0, zaczynając od wartości przybliżonej (bez oporu powietrza) *)
X=x/.rozwiazanierr;
Y=y/.rozwiazanierr;
{X[tZ],Y[tZ]} (* po wstawieniu czasu lotu otrzymujemy poszukiwane współrzędne upadku pocisku *)
]

wspolrzedne[1,3,1,1,10,1,0.1] (* test *)

(* Funkcja zwracająca losowy wektor prędkości o długości mniejszej niż Vmax, patrz Zad 4, Zestaw 6 *)
predkoscLosowa[Vmax_]:=Module[{n},
Label[poczatek];
n=RandomReal[{-1,1},3];
If[n.n≤1 ,Return[n Vmax],Goto[poczatek]]
]

Table[predkoscLosowa[0.8],{10000}];

ListPointPlot3D[%,PlotRange->{{-1,1},{-1,1},{-1,1}},BoxRatios->{1,1,1},PlotStyle->PointSize[Tiny],RotationAction->"Clip"] (* test losowania prędkości *)

cel[{x_,y_}]:=If[Abs[x-3]≤1&&Abs[y-4]≤2,True,False]
(* przykładowy cel: prostokąt *)

cel[{2.4,2.3}]

RegionPlot[cel[{x,y}],{x,0,10},{y,0,10}] (* celujemy w ten kwadrat *)

(* Własciwa symulacja Monte Carlo *)

symulacja1={};
symulacja2={};
For[k=1,k≤2000,k++,
{V0x,V0y,V0z}=predkoscLosowa[7];
hit=wspolrzedne[V0x,V0y,V0z,10,9.81,1,0.2];
If[cel[hit],AppendTo[symulacja1,hit],AppendTo[symulacja2,hit]]
]

trafione=Dimensions[symulacja1]//First

nietrafione=Dimensions[symulacja2]//First

Show[ListPlot[symulacja1,PlotStyle->Red],ListPlot[symulacja2,PlotStyle->Blue],PlotRange->All,AspectRatio->1,AxesOrigin->{0,0}](* wizualizacja trafień *)

prawdopodobienstwo=trafione/(trafione+nietrafione)//N