Post Reply 
Earthquake: Earthquake Calculations
07-15-2019, 10:54 AM (This post was last modified: 07-15-2019 01:05 PM by StephenG1CMZ.)
Post: #6
RE: Earthquake: Earthquake Calculations
Version 1.2:
Estimates M0 and magnitude.
Estimates distance/time.
Compares relative size and strength.
Most functions allow lists of values.
In English, francaise and polskie.

Code:
 
 
 LOCAL VER:=" V1.2: ";
 LOCAL MN:="Earthquake"+VER; 
 LOCAL CRID:=MN+" © 2019 StephenG1CMZ";
 
 LOCAL LANGSE:={" English ",""," French ","","","","",""," Polish "};
 LOCAL LANGS :={" English ",""," français ","","","","",""," polskie "};
 LOCAL LANGS2:={"EN","","FR","","","","","","PL"}; //EN CHINESE FR GERMAN SPANISH DUTCH PORTUGUESE JAPANESE PL
 LOCAL SL:=1;//SL

 EXPORT RigidityTypical:=3ᴇ10; //FOR CONVENIENCE
 
 //VELOCITY in unit/s (to match Centre distance)
 EXPORT VP; //TYPICAL 6 km/s OR 1.7*VS
 EXPORT VS; //TYPICAL 3 km/s
   //THE TYPICAL VALUES HERE REFER TO LOCAL WAVES
  
 LOCAL EPSILON:=1ᴇ−308;//STD MATH:SMALLEST NONZERO (LIMIT VALUE USED FOR LOG0)
 LOCAL MNF:={"Earthquake","","Tremblement de terre","","","","","","Trzęsienie ziemi"}+VER;
 LOCAL MAGT:={"Magnitude","","Grandeur","","","","","","Wielkość"};
 LOCAL RELT:={"Relative ","","Relative ","","","","","","względna "};
 LOCAL AMPLT:={"Amplitude","","Amplitude","","","","","","Amplituda"};
 LOCAL ENERGT:={"Energy","","Énergie","","","","","","Energia"};
 ABOUT()
 BEGIN
  
  LOCAL ABT:={"Earthquake calculations are estimates and can be calculated differently. No liability is accepted.","",
   "Les calculs sismiques sont des estimations et peuvent être calculés différemment.  Aucune responsabilité n'est acceptée.","","","","","",
   "Obliczenia trzęsienia ziemi są szacunkowe i można je obliczyć inaczej.  Nie ponosimy odpowiedzialności."};
  
  PRINT();
  PRINT(CRID);
  PRINT(MNF(SL));
  PRINT(ABT(SL));
  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 lists
 
 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:=IFTE(RigidityNm2,RigidityNm2,RigidityTypical);//Use Typical value if none supplied
  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 ROUND(ALOG(5.24+1.44*Magnitud),0);//J
 END;

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

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

 //Timing requires VP and VS to be set
 //Typically km/s but ensure these units match
 //VP VS (km/s) DISTANCE (km) TDELTA (s)
 
 EXPORT GetVPVS()
 BEGIN
 //PROMPT USER FOR VP VS 
 //(CALLABLE BY USER OR AUTOMAGICALLY TO AVOID ERROR)
 //TYPICALLY 6km/s and ?
  LOCAL TTL:={"Enter Velocity of P and S waves","","Entrez la vitesse des ondes P et S","","","","","","Wprowadź prędkość fal P i S"};
  LOCAL LBL:={"VP","VS"};
  LOCAL HLPP:={"Primary wave (typically  6 km/s, 1.7*VS)","","Vague primaire (généralement 6 km / s, 1.7*VS)","","","","","","Fala pierwotna (zazwyczaj 6 km / s, 1.7*VS)"};
  LOCAL HLPS:={"Secondary wave (typically  3 km/s)","","vague secondaire (généralement 6 km / s)","","","","","","Fala wtórna (zazwyczaj 3 km / s)"}; 
  LOCAL HLP:={HLPP(SL),HLPS(SL)};
  INPUT({VP,VS},TTL(SL),LBL,HLP);
 END;

 EXPORT CentreL(TDELTA)
 //Distance to centre-Local usable up To 100-250km.
 //Simple formula for constant velocity
 //TDELTA=(TimeS-TimeP) s (alt. to match VP VS)
 //RESULTS in m or km (ie. as VP VS) 
 BEGIN
  LOCAL DST;
  IF VP-VS==0 OR VP==0 OR VS==0 THEN
   MSGBOX(MN+"CentreL\n"+"VP-VS=0");
   GetVPVS();
  END;
  //MSGBOX({TDELTA,VP,VS});//INPUTLOG
  DST:=(TDELTA)*(VP*VS)/(VP-VS);
  RETURN ROUND(DST,0);//unit matching VS//ROUND to km
 END;
 //These trivial distace-time-velocity formulae assume constant velocity 
 EXPORT Delay(DISTANCE) // DISTANCE ≤250 km
 //GIVEN DISTANCE,HOW LONG AFTER P WILL S ARRIVE
 //DISTANCE:UNITS TO MATCH VP VS 
 //RESULTS MATCH UNITS OF VELOCITIES VP VS
 BEGIN
  IF VP==0 OR VS==0 THEN
   MSGBOX(MN+"\n"+"VP/0 VS/0:"+STRING({VP,VS}));
   GetVPVS();
  END;
  //MSGBOX({DISTANCE,VP,VS});//INPUTLOG
  RETURN (DISTANCE/VP)-(DISTANCE/VS);
 END;
 
 EXPORT AnEarthquake(Magnitud)
 //Describe an Earthquake
 BEGIN
  PRINT();
  PRINT(STRING(Magnitud)+" Magnitude");
  PRINT(" ");
  PRINT(STRING(EnergyJ(Magnitud))+" J");
  PRINT(STRING(Energyt(Magnitud))+" t(TNT)");  
  PRINT(STRING(EnergyH(Magnitud))+" Hiroshima");
 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 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;

 EXPORT Earthquakes(Magnitudes,INDX)
 //Describe relative earthquakes
 //Values onscreen need to be rounded to 9
 //to prevent occasional overlap. 
 //You may prefer to omit this rounding.
 BEGIN
  LOCAL II;
  LOCAL DP:=9;
  LOCAL SW:=320;
  LOCAL AMPLS:=Amplitudes(Magnitudes,INDX);
  LOCAL STRONGS:=Energies(Magnitudes,INDX);

  RECT_P();
  TEXTOUT_P(MN+MNF(SL),0,220); 
  TEXTOUT_P(MAGT(SL),0,0);
  TEXTOUT_P(AMPLT(SL),SW/5,0);
  TEXTOUT_P(ENERGT(SL),2*SW/4,0);
  TEXTOUT_P(RELT(SL)+INDX,3*SW/4,0); 
  IF TYPE(Magnitudes)==TYPE({}) THEN
   FOR II FROM 1 TO MIN(SIZE(Magnitudes),10) DO //10-11 PER SCREEN:1SCREEN 
    TEXTOUT_P(ROUND(Magnitudes(II),4),0,20*II); 
    TEXTOUT_P(ROUND(AMPLS(II),DP),SW/6,20*II);     //ROUND FOR DISPLAY
    TEXTOUT_P(ROUND(STRONGS(II),DP),2*SW/4,20*II); //ROUND FOR DISPLAY
   END;
  END;
  FREEZE;
  WAIT;
 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 ");
  WAIT;
  AnEarthquake(7.1);
 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");
  WAIT;
  Earthquakes({1,2,3,4,5,6,7,8,9,10,11,12,13},4);
 END;

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

 EXPORT HELP()
 BEGIN
  LOCAL ALL:={"Most functions allow lists","","La plupart des fonctions permettent les listes","","","","","","Większość funkcji pozwala na listy"};
  LOCAL MAG:=MAGT(SL)+": "+{"Not exactly","","Pas exactement","","","","","","Nie dokładnie"}+" Richter";
  LOCAL UL:= {"Unitless","","Sans unité","","","","","","Bez jednostki"};
  ABOUT();
  PRINT();
  PRINT(ALL(SL));
  PRINT(" ");
  //PRINT("UNITS");
  PRINT("M0Nm: Nm");
  PRINT(MAGT(1)+": "+MAG(SL));
  PRINT(ENERGT(SL)+": J, Hiroshimas or tonnesTNT");
  PRINT("GetVPVS: km/s?");
  PRINT("CentreL: km  ? (≤100-250km");
  PRINT("Delay:      s?");
  PRINT(" ");
  PRINT("RelativeAmplitudes: "+UL(SL));
  PRINT("RelativeEnergies: "+UL(SL));
  PRINT(" ");
 END;

 EXPORT Languages()
 BEGIN
  LOCAL TTL:=IFTE(SL,LANGS2(SL)+LANGS(SL)+LANGSE(SL),LANGS2+LANGSE); 
  IF CHOOSE(SL,TTL,LANGS2+LANGS+LANGSE) THEN
   MSGBOX(LANGS2(SL)+LANGS(SL)+LANGSE(SL));
  END;
 END;

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


References providing formulae and examples
1. https://www.e-education.psu.edu/earth520...l7_p4.html - Earthquake magnitude and energy
2. https://www.hpmuseum.org/forum/thread-13239.html - Relative strength of 2 earthquakes
3. https://www.google.com/url?sa=t&source=w...NFcE6pGX2B" - PDF - Earthquake location and speed

Stephen Lewkowicz (G1CMZ)
https://my.numworks.com/python/steveg1cmz
Visit this user's website Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
RE: Earthquake: Earthquake Calculations - StephenG1CMZ - 07-15-2019 10:54 AM



User(s) browsing this thread: 1 Guest(s)