Advertisement
NUCLEARESOL

Finding least time for curves

Dec 3rd, 2023
814
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C# 16.31 KB | None | 0 0
  1. using System;
  2. using System.Diagnostics;
  3. using System.Threading;
  4. using System.Globalization;
  5.  
  6. class Program
  7. {
  8.  
  9.     static void Main()
  10.     {
  11.         int version = 1;
  12.         string type = "circle";
  13.         double gravitationalA = 9.81d;
  14.         double[] startPoint = { 5,10};
  15.         double[] endpoint = { 0, 0 };
  16.         const Boolean calculateDistance = false;
  17.  
  18.         Console.WriteLine("Following Variables are used:");
  19.         Console.WriteLine("Gravitational Accel: " + gravitationalA + " m/s");
  20.         Console.WriteLine("Goal: " + endpoint[0] + "," + endpoint[1]);
  21.         Console.WriteLine("Start: " + startPoint[0] + "," + startPoint[1]);
  22.         Console.WriteLine("Calculate Distance" + calculateDistance);
  23.         Console.WriteLine("y/n");
  24.         string confirmation = Console.ReadLine();
  25.  
  26.         if (confirmation != "y")
  27.         {
  28.             Console.WriteLine("");
  29.             Console.WriteLine("Abort, Please change the variable");
  30.             Console.WriteLine("graph type");
  31.             type = Console.ReadLine();
  32.             Console.WriteLine("Gravitational Acceleration");
  33.             gravitationalA = double.Parse(Console.ReadLine(), CultureInfo.InvariantCulture.NumberFormat);
  34.             Console.WriteLine("Start Point X");
  35.             startPoint[0] = double.Parse(Console.ReadLine(), CultureInfo.InvariantCulture.NumberFormat);
  36.             Console.WriteLine("Start Point Y");
  37.             startPoint[1] = double.Parse(Console.ReadLine(), CultureInfo.InvariantCulture.NumberFormat);
  38.             Console.WriteLine("");
  39.             Console.WriteLine("Following Variables are used:");
  40.             Console.WriteLine("Gravitational Accel: " + gravitationalA + " m/s");
  41.             Console.WriteLine("Goal: " + endpoint[0] + "," + endpoint[1]);
  42.             Console.WriteLine("Start: " + startPoint[0] + "," + startPoint[1]);
  43.             //Thread.Sleep(3000);
  44.         }
  45.  
  46.         // Console.WriteLine("");
  47.         // Console.WriteLine("Starting calculation");
  48.         double x = startPoint[0];
  49.         double velX = 0f;
  50.         const double timestep = 0.0001f;
  51.         double time = 0f;
  52.         double initialCollvar = 1;
  53.         double initialCollvarStep = 0.1f;
  54.         double zeroInAccuracy = 0.0001f;
  55.         double maxTime = double.PositiveInfinity;
  56.         double limit = 400;
  57.         string strg = "";
  58.         Console.WriteLine(limit);
  59.         var start = Process.GetCurrentProcess().TotalProcessorTime;
  60.         int i = 0;
  61.         static double NewtonsMethod(double target, double initialGuess, double tolerance, int maxIterations)
  62.         {
  63.             double x = initialGuess;
  64.  
  65.             for (int i = 0; i < maxIterations; i++)
  66.             {
  67.                 double fx = (2 * x - Math.Sin(2 * x)) / (1 - Math.Cos(2 * x));
  68.                 double derivative = -2 * ((x / Math.Tan(x)) - 1) * Math.Pow(1 / Math.Sin(x), 2);
  69.                 Console.WriteLine(x + " " + fx + " " + derivative);
  70.                 //Thread.Sleep(1000);
  71.                 if (Math.Abs(fx - target) < tolerance)
  72.                 {
  73.                     Console.WriteLine($"Converged after {i} iterations.");
  74.                     return x;
  75.                 }
  76.  
  77.                 x = x - (fx - target) / derivative;
  78.             }
  79.  
  80.             Console.WriteLine("Did not converge within the specified number of iterations.");
  81.             return x;
  82.         }
  83.  
  84.         Tuple<double, double> calculateAcceleration(double X, double collVar, double depvar1, double depvar2)
  85.         {
  86.             double localSlope = 2.0f * collVar * (X - depvar1);
  87.  
  88.             //Console.WriteLine(-(float)(Math.Cos(Math.Atan(1f / localSlope)) * Math.Sin(Math.Atan(1f / localSlope)) * gravitationalA));
  89.             return Tuple.Create(-(double)(Math.Cos(Math.Atan(1f / localSlope)) * Math.Sin(Math.Atan(1f / localSlope)) * gravitationalA), localSlope);
  90.         }
  91.         Tuple<double, double> calculateAccelerationCircle(double X, double Y, double collVar, double h, double k)
  92.         {
  93.             double localSlope = (h - X) / (Y - k);
  94.  
  95.             //return -(float)(Math.Cos(Math.Atan(1f / localSlope)) * Math.Sin(Math.Atan(1f / localSlope)) * gravitationalA);
  96.             return Tuple.Create(-(double)(Math.Cos(Math.Atan(1f / localSlope)) * Math.Sin(Math.Atan(1f / localSlope)) * gravitationalA), localSlope);
  97.         }
  98.         Tuple<double, double> calculateAccelerationCubic(double X, double collVar, double h, double k)
  99.         {
  100.             double localSlope = collVar*3.0d*Math.Pow(X-h,2.0);
  101.  
  102.             //return -(float)(Math.Cos(Math.Atan(1f / localSlope)) * Math.Sin(Math.Atan(1f / localSlope)) * gravitationalA);
  103.             return Tuple.Create(-(double)(Math.Cos(Math.Atan(1f / localSlope)) * Math.Sin(Math.Atan(1f / localSlope)) * gravitationalA), localSlope);
  104.         }
  105.  
  106.         Console.WriteLine(strg);
  107.         Console.WriteLine("");
  108.         strg = "";
  109.         Tuple<double, double, double, double> CalculateTimeParabola(double initialCollvar)
  110.         {
  111.             double h = 0;
  112.             double k = 0;
  113.             double distance = 0;
  114.             double deltaX = 0;
  115.             i = i + 1;
  116.             x = startPoint[0];
  117.             velX = 0f;
  118.             time = 0f;
  119.             distance = 0f;
  120.             h = 0.5f * (startPoint[0] - (startPoint[1] / (startPoint[0] * initialCollvar)));
  121.             k = -initialCollvar * (double)Math.Pow(h, 2);
  122.             while (x >= 0 && time < limit)
  123.             {
  124.                 Tuple<double, double> calculatedVales = calculateAcceleration(x, initialCollvar, h, k);
  125.                 deltaX = (double)(velX * timestep + 0.5f * calculatedVales.Item1 * Math.Pow(timestep, 2));
  126.                 x = (double)(x + deltaX);
  127.                 velX = velX + (0 + calculateAcceleration(x, initialCollvar, h, k).Item1) * timestep / 1f;
  128.                 if (calculateDistance)
  129.                 {
  130.                     distance = distance + (double)Math.Sqrt(deltaX * deltaX * (calculatedVales.Item2 * calculatedVales.Item2 + 1.0f));
  131.                 }
  132.                 time = time + timestep;
  133.             }
  134.             Console.WriteLine(Math.Abs(Math.Abs(initialCollvarStep) * 100 - 100) + " % complete");
  135.  
  136.             return Tuple.Create(time, initialCollvar, h, k);
  137.         }
  138.  
  139.         Tuple<double, double, double, double> CalculateTimeCircle(double initialCollvar)
  140.         {
  141.             if (initialCollvar < (double)Math.Sqrt(Math.Pow(startPoint[1], 2) + Math.Pow((Math.Pow(startPoint[1], 2) - Math.Pow(startPoint[0], 2)) / (-2 * startPoint[0]), 2)))
  142.             {
  143.                 return Tuple.Create(limit, initialCollvar, 0.0d, 0.0d);
  144.             }
  145.             double q = Math.Sqrt(Math.Pow(startPoint[1], 2) + Math.Pow(startPoint[0], 2));
  146.             double x3 = startPoint[0] / 2;
  147.             double y3 = startPoint[1] / 2;
  148.  
  149.  
  150.             double h = x3 + (Math.Sqrt(Math.Pow(initialCollvar, 2) - Math.Pow(q / 2, 2)) * (-startPoint[1])) / q;
  151.             double k = y3 + (Math.Sqrt(Math.Pow(initialCollvar, 2) - Math.Pow(q / 2, 2)) * (startPoint[0])) / q;
  152.  
  153.             double distance = 0;
  154.             double deltaX = 0;
  155.  
  156.             x = startPoint[0];
  157.             velX = 0f;
  158.             time = 0f;
  159.             distance = 0f;
  160.             double y = k - Math.Sqrt(Math.Pow(initialCollvar, 2) - Math.Pow(x - h, 2));
  161.             while (x >= 0 && time < limit)
  162.             {
  163.                 y = k - Math.Sqrt(Math.Pow(initialCollvar, 2) - Math.Pow(x - h, 2));
  164.                 Tuple<double, double> calculatedVales = calculateAccelerationCircle(x, y, initialCollvar, h, k);
  165.                 deltaX = (double)(velX * timestep + 0.5f * calculatedVales.Item1 * Math.Pow(timestep, 2));
  166.                 x = (double)(x + deltaX);
  167.                 velX = velX + (0 + calculateAccelerationCircle(x, y, initialCollvar, h, k).Item1) * timestep / 1f;
  168.                 if (calculateDistance)
  169.                 {
  170.                     distance = distance + (double)Math.Sqrt(deltaX * deltaX * (calculatedVales.Item2 * calculatedVales.Item2 + 1.0f));
  171.                 }
  172.                 time = time + timestep;
  173.             }
  174.  
  175.             return Tuple.Create(time, initialCollvar, h, k);
  176.         }
  177.         Tuple<double, double, double, double> CalculateTimeCubic(double initialCollvar)
  178.         {
  179.             double h = 0;
  180.             double k = 0;
  181.             double distance = 0;
  182.             double deltaX = 0;
  183.             i = i + 1;
  184.             x = startPoint[0];
  185.             velX = 0f;
  186.             time = 0f;
  187.             distance = 0f;
  188.             if (initialCollvar <= 0.00001 |initialCollvar >= 12.0 * startPoint[1] / (3.0 * Math.Pow(startPoint[0], 3)))
  189.             {
  190.                 return Tuple.Create(limit+100, initialCollvar, h, k);
  191.             }
  192.  
  193.             h = (3.0d * Math.Pow(startPoint[0], 2) - Math.Sqrt(9.0d * Math.Pow(startPoint[0], 4) - 4.0 * 3.0 * startPoint[0] * (Math.Pow(startPoint[0], 3) - startPoint[1] / initialCollvar)))/2;
  194.             k = initialCollvar * (double)Math.Pow(h, 3);
  195.             while (x >= 0 && time < limit)
  196.             {
  197.                 Tuple<double, double> calculatedVales = calculateAccelerationCubic(x, initialCollvar, h, k);
  198.                 deltaX = (double)(velX * timestep + 0.5f * calculatedVales.Item1 * Math.Pow(timestep, 2));
  199.                 x = (double)(x + deltaX);
  200.                 velX = velX + (0 + calculateAccelerationCubic(x, initialCollvar, h, k).Item1) * timestep / 1f;
  201.                 if (calculateDistance)
  202.                 {
  203.                     distance = distance + (double)Math.Sqrt(deltaX * deltaX * (calculatedVales.Item2 * calculatedVales.Item2 + 1.0f));
  204.                 }
  205.                 time = time + timestep;
  206.             }
  207.             Console.WriteLine(Math.Abs(Math.Abs(initialCollvarStep) * 100 - 100) + " % complete");
  208.  
  209.             return Tuple.Create(time, initialCollvar, h, k);
  210.         }
  211.  
  212.         Tuple<double, double, double, double, double, double> jumpingDecend()
  213.         {
  214.             //   Console.WriteLine(Math.Abs(initialCollvarStep)+","+zeroInAccuracy);
  215.             double var1 = 0;
  216.             double var2 = 0;
  217.             double var3 = 0;
  218.             double var4 = 0;
  219.             if (type == "parabola")
  220.             {
  221.                 var1 = 0;
  222.                 var2 = 0;
  223.             }
  224.             if (type == "circle")
  225.             {
  226.                 initialCollvar = (double)Math.Sqrt(Math.Pow(startPoint[1], 2) + Math.Pow((Math.Pow(startPoint[1], 2) - Math.Pow(startPoint[0], 2)) / (-2 * startPoint[0]), 2));
  227.                 //  Console.WriteLine(initialCollvar+"RADIUS");
  228.             }
  229.             if (type == "cubic")
  230.             {
  231.                 initialCollvar = (12.0 * startPoint[1] / (3.0 * Math.Pow(startPoint[0], 3)))/2;
  232.                 initialCollvarStep = (12.0 * startPoint[1] / (3.0 * Math.Pow(startPoint[0], 3))) / 10;
  233.                 zeroInAccuracy = (12.0 * startPoint[1] / (3.0 * Math.Pow(startPoint[0], 3))) / 100000;
  234.                 Console.WriteLine(initialCollvar+"RADIUS");
  235.             }
  236.  
  237.             while (Math.Abs(initialCollvarStep) >= zeroInAccuracy)
  238.             {
  239.                 // --  Console.Write("BOB");
  240.                 initialCollvar = initialCollvar + initialCollvarStep;
  241.                 if (type == "parabola")
  242.                 {
  243.                     time = CalculateTimeParabola(initialCollvar).Item1;
  244.                 }
  245.                 else if (type == "circle")
  246.                 {
  247.                     time = CalculateTimeCircle(initialCollvar).Item1;
  248.                 }
  249.                 else if (type == "cubic")
  250.                 {
  251.                     time = CalculateTimeCubic(initialCollvar).Item1;
  252.                 }
  253.                 // Console.WriteLine(distance);
  254.                 //  Console.WriteLine(initialCollvar + "(x-" + h + ")^2+" + k);
  255.                 strg = strg + "(" + initialCollvar + "," + time + ")" + ",";
  256.                 if (maxTime > time)
  257.                 {
  258.                     maxTime = time;
  259.                     initialCollvarStep = initialCollvarStep * 2;
  260.                 }
  261.  
  262.                 else if (time > maxTime)
  263.                 {
  264.                     if (time >= limit)
  265.                     {
  266.                         initialCollvar = initialCollvar - initialCollvarStep;
  267.                         initialCollvarStep = -initialCollvarStep / 4;
  268.                     }
  269.                     else
  270.                     {
  271.                         initialCollvarStep = -initialCollvarStep;
  272.                     }
  273.                 }
  274.                 else
  275.                 {
  276.                     initialCollvarStep = initialCollvarStep / 2;
  277.                 }
  278.             }
  279.             if (type == "circle")
  280.             {
  281.                 Tuple<double, double, double, double> resultcal = CalculateTimeCircle(initialCollvar);
  282.                 time = CalculateTimeCircle(initialCollvar).Item1;
  283.                 var1 = resultcal.Item1;
  284.                 var2 = resultcal.Item3;
  285.                 var3 = resultcal.Item4;
  286.             }
  287.             else if (type == "parabola")
  288.             {
  289.                 Tuple<double, double, double, double> resultcal = CalculateTimeParabola(initialCollvar);
  290.                 time = CalculateTimeParabola(initialCollvar).Item1;
  291.                 var1 = resultcal.Item2;
  292.                 var2 = resultcal.Item3;
  293.                 var3 = resultcal.Item4;
  294.             }
  295.             else if (type == "cubic")
  296.             {
  297.                 Tuple<double, double, double, double> resultcal = CalculateTimeCubic(initialCollvar);
  298.                 time = CalculateTimeCubic(initialCollvar).Item1;
  299.                 var1 = resultcal.Item2;
  300.                 var2 = resultcal.Item3;
  301.                 var3 = resultcal.Item4;
  302.             }
  303.             if (time > limit)
  304.             {
  305.  
  306.                 if (type == "circle")
  307.                 {
  308.                     initialCollvar = initialCollvar - initialCollvarStep / 2;
  309.                     Tuple<double, double, double, double> resultcal = CalculateTimeCircle(initialCollvar);
  310.                     time = resultcal.Item1;
  311.                     var1 = resultcal.Item2;
  312.                     var2 = resultcal.Item3;
  313.                     var3 = resultcal.Item4;
  314.                 }
  315.                 else if (type == "parabola")
  316.                 {
  317.                     initialCollvar = initialCollvar - initialCollvarStep / 2;
  318.                     Tuple<double, double, double, double> resultcal = CalculateTimeParabola(initialCollvar);
  319.                     time = resultcal.Item1;
  320.                     var1 = resultcal.Item2;
  321.                     var2 = resultcal.Item3;
  322.                     var3 = resultcal.Item4;
  323.                 }
  324.                 else if (type == "cubic")
  325.                 {
  326.                     initialCollvar = initialCollvar - initialCollvarStep / 2;
  327.                     Tuple<double, double, double, double> resultcal = CalculateTimeCubic(initialCollvar);
  328.                     time = resultcal.Item1;
  329.                     var1 = resultcal.Item2;
  330.                     var2 = resultcal.Item3;
  331.                     var3 = resultcal.Item4;
  332.                 }
  333.             }
  334.             Console.WriteLine("Finished");
  335.             return Tuple.Create(initialCollvar, time, var1, var2, var3, var4);
  336.         }
  337.         Tuple<double, double, double, double, double, double> resultTuple = jumpingDecend();
  338.         var stop = Process.GetCurrentProcess().TotalProcessorTime;
  339.  
  340.         Console.WriteLine("");
  341.         Console.WriteLine("time taken " + resultTuple.Item2 + " seconds");
  342.         Console.WriteLine("At variable = " + resultTuple.Item1);
  343.         Console.WriteLine("h = " + resultTuple.Item4);
  344.         Console.WriteLine("k = " + resultTuple.Item5);
  345.         Console.WriteLine();
  346.         // stopwatch.Stop();
  347.         Console.WriteLine("");
  348.         Console.WriteLine("Calculation Time is " + (stop - start).TotalMilliseconds);
  349.         Console.WriteLine("Cal amt" + i);
  350.       //  Console.WriteLine("Starting control group in 5 seconds");
  351.         Console.WriteLine(strg);
  352.         Thread.Sleep(20000);
  353.         start = Process.GetCurrentProcess().TotalProcessorTime;
  354.         double theta = NewtonsMethod(startPoint[0] / startPoint[1],1,1e-20,2000);
  355.         double D = 2 * startPoint[1] / (1 - Math.Cos(theta * 2));
  356.         Console.WriteLine(D);
  357.         double t = Math.Sqrt(2*D/gravitationalA)*theta;
  358.         Console.WriteLine(t);
  359.         stop = Process.GetCurrentProcess().TotalProcessorTime;
  360.         Console.WriteLine("Calculation Time is " + (stop - start).TotalMilliseconds);
  361.         strg = "";
  362.     }
  363. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement