Earthquake: Earthquake Calculations

[StephenG1CMZ]

Earthquake V 1.1
All functions can now handle lists of values.
Error handling improved.
A list of values can be compared for relative amplitude and energy.

Code:


 LOCAL MN:="Earthquake V1.1: ";
 LOCAL CRID:=MN+" © 2019 StephenG1CMZ"; 

 EXPORT RigidityTypical:=3ᴇ10; //FOR CONVENIENCE

 LOCAL EPSILON:=1ᴇ−308;//SMALLEST NONZERO (LIMIT VALUE USED FOR LOG0)
 
 ABOUT()
 BEGIN
  PRINT();
  PRINT(CRID);
  PRINT("Earthquake calculations are estimates");
  PRINT("and can be calculated differently.");
  PRINT("No liabilIty is accepted");
  WAIT;
 END;

 //STD MATH
 LOG10(NUM) //AVOID RISK OF MISTAKING LOGS
 BEGIN
  RETURN LOG(NUM);//LOG10 
 END;
 //END MATH

 // Calculations for a single quake
 // Parameters may be numeric or list
 
 EXPORT AnEarthquake__() //HEADER
 BEGIN

 END;

 EXPORT M0Nm(RigidityNm2,Aream2,Slipm)
 //Calculate SeismicMoment in Nm^2
 //=rigidity x area x slip
 //The rigidity of rock is a constant number based on the rock type. It has units of pressure. Typical assumptions are on the order of 3 x 10^10 N/m2.
 //Slip is a length and it is on the order of centimeters (meters for a great earthquake). 
 //Area is in units of length2 and is often on the order of km2. 
 //The units for seismic moment are then Nm (newton meters). 
 //As an example, the 2004 26 December Sumatra-Andaman earthquake had the following dimensions as reported by Lay et al. (2005): Its slip averaged about 5 m, its rupture length was about 1300 km and the fault width was between 160 - 240 km. Assuming a rigidity of 3 x 10^10 N/m2 gives us a seismic moment of 4.4 x 10^22Nm.
 BEGIN
  LOCAL Rigidity:=RigidityNm2;
  LOCAL Area:=Aream2;
  LOCAL Slip:=Slipm;
  RETURN Rigidity*Area*Slip;
 END;

 MW(M0)
 BEGIN
  //Calculate Mw = (2/3)*logM0 - 6.05
  LOCAL M00:=MAX(M0,EPSILON);//GUARD M0≤0
  LOCAL MM:=(2/3)*LOG10(M00)-6.05;
  IF MIN(M0)≤0 THEN
   MSGBOX(MN+"MW:\n"+"LOG≤0: "+STRING(M0));
  END; 
  RETURN MM;
 END;

 EXPORT Magnitude(RigidityNm2,Aream2,Slipm)
 BEGIN
  RETURN MW(M0Nm(RigidityNm2,Aream2,Slipm));
 END;
 //These magnitudes are approximations to Richter scale but not same
 EXPORT Magnitude_(M0)
 BEGIN
  RETURN MW(M0);
 END;

 EXPORT EnergyJ(Magnitud) // //"Magnitude" is a reserved word or function
 // Calculate E, where Log10(E) = 5.24 + 1.44M
 //Empirical good for M>5
 //This estimates damaging energy
 BEGIN
  RETURN ALOG(5.24+1.44*Magnitud);//J
 END;

 EXPORT EnergyH(Magnitud)
 BEGIN
  LOCAL HIROSHIMA:=7.4ᴇ12;//J
  RETURN EnergyJ(Magnitud)/HIROSHIMA;
 END;

 EXPORT Energyt(Magnitud)
 BEGIN
  LOCAL TNT:=4.184ᴇ9;//J PER TONNE
  RETURN EnergyJ(Magnitud)/TNT;
 END;

 //Comparing 2 Quakes (unused)

 RelAmplitude(MagA,MagB)
 //Compare relative amplitude, any order 
 BEGIN
  LOCAL BIGGER := 10^ABS(MagA-MagB);
  RETURN BIGGER; //unitless
 END;

 RelEnergy(MagA,MagB)
 //Compare relative energy (damage), any order
 BEGIN
  LOCAL BIGGER   := RelAmplitude(MagA,MagB);
  LOCAL STRONGER := 10^(ABS(MagA-MagB)/2)*BIGGER;
  RETURN STRONGER; //unitless
 END;
 
 // Comparing a list of quakes

 EXPORT ComparingEarthquakes__() //HEADER
 BEGIN

 END;

 EXPORT Amplitudes(Magnitudes,INDX)
 //Compare relative amplitude with indexed value=1
 //Errors return empty list (except LOG:LIMITED)
 BEGIN
  LOCAL MAG;
  LOCAL BIGGER;
  
  CASE
   IF TYPE(Magnitudes)=Type({}) THEN
    IF SIZE(Magnitudes)=0 OR INDX>SIZE(Magnitudes) THEN
     RETURN {};
    END;
    MAG:=Magnitudes(INDX);
    BIGGER := 10^(Magnitudes-MAG);
    RETURN BIGGER; //unitless
   END;  
   DEFAULT
    RETURN {}; //{MN+"Provide a list"};
  END; 
 END;

 EXPORT Energies(Magnitudes,INDX)
 //Compare relative energy (damage) with indexed value=1
 //Errors return empty list (except LOG:LIMITED)
 BEGIN
  LOCAL MAG;
  LOCAL BIGGER;
  LOCAL STRONGER;
   CASE
    IF TYPE(Magnitudes)=Type({}) THEN
     IF SIZE(Magnitudes)=0 OR INDX>SIZE(Magnitudes) THEN
      RETURN {};
     END;
     MAG:=Magnitudes(INDX);
     BIGGER   := Amplitudes(Magnitudes,INDX);
     STRONGER := 10^((Magnitudes-MAG)/2)*BIGGER; 
     RETURN STRONGER; //unitless
    END;
    DEFAULT  
     RETURN {}; //{MN+"Provide a list"};
   END;
 END;
 
 //Examples

 ExM0()
 BEGIN
  // Calculate upper and lower bounds for ex
  // Demonstrates mixing list and numeric
  PRINT ({M0Nm(3*10^10,{1300ᴇ3*(160ᴇ3),1300ᴇ3*(240ᴇ3)},5)," =4.4 x 10^22"}+" Nm ");
 END;

 ExMW()
 BEGIN
  PRINT ({MW(4.4ᴇ22),"=9.05"}+" Magnitude ");
 END;

 ExEnergy()
 BEGIN
  PRINT({EnergyJ({5.0,6.0,7.0}),"=2.8ᴇ12 =7.8ᴇ13?? =2.1ᴇ15"}+" J ");
  PRINT({EnergyH(5.0),EnergyH(6.0),EnergyH(7.0)," "}+" Hiroshima ");
  PRINT({EnergyJ(9.05),"=1.84ᴇ18"}+" J ");
  PRINT({Energyt(5.0)}+" tTNT ");
 END;

 ExRel()
 BEGIN
  PRINT({RelAmplitude(7.1,6.4)," =5"});
  PRINT({RelEnergy   (7.1,6.4)," =11.2"});
  //PRINT(Amplitudes({7.1,6.4},1));
  //PRINT(Energies({7.1,6.4},1));
  PRINT(STRING(Amplitudes({7.1,6.4},2))+" =5");
  PRINT(STRING(Energies({7.1,6.4},2))+" =11.2");
 END;

 EXPORT Examples()
 BEGIN 
 
  ABOUT();
  //PRINT();
  PRINT("[=] signifies expected value");
  WAIT; 
  ExM0();
  ExMW();
  ExEnergy();
  ExRel();
 END;

 EXPORT HELP()
 BEGIN
  ABOUT();
  PRINT();
  PRINT("All functions allow lists");
  PRINT("UNITS");
  PRINT("M0Nm: M0 in Nm");
  PRINT("Magnitude: not exactly Richter");
  PRINT("Energy: in J, Hiroshimas or tonnesTNT");
  PRINT("RelativeAmplitudes: unitless");
  PRINT("RelativeEnergies: unitless ");
  PRINT(" ");
 END;

 EXPORT References()
 BEGIN
  PRINT("formulae and examples-
https://www.e-education.psu.edu/earth520/content/l7_p4.html");
  PRINT("https://www.hpmuseum.org/forum/thread-13239.html");
 END;

 EXPORT EARTHQUAKE() //MAIN
 BEGIN
   ABOUT();
 END;