// Syntax to run from Stata: framingham30, male(male) age(age) sbp(sbp) smoke(smoke) trtbp(trtbp) diab(diab) bmi(bmi) tcl(tcl) hdl(hdl)

//---------------------------- Main Stata program --------------------------

program drop _all

program define framingham30, eclass
version 13
	syntax [if] [in] ,  ///
		Male(varname numeric)  	/// integer 0 (female) or 1 (male)
		AGE(varname numeric)    /// age in years
		SBP(varname numeric)  	/// Systolic Blood Pressure (SBP)
		SMoke(varname numeric) 	/// integer 0 (non-smoking) or 1 (smoking)
		TRtbp(varname numeric) 	/// Use of Antihypertensive treatment: integer 0 (no) or 1 (yes)
		DIAB(varname numeric) 	/// Diabetes mellitus: integer 0 (no) or 1 (yes)
		[BMI(varname numeric) REPLace SUFFix(str) *] 	/// Body mass index
		[TCL(varname numeric) REPLace SUFFix(str) *]  	/// Total cholesterol
		[HDL(varname numeric) REPLace SUFFix(str) *]  	/// HDL cholesterol

marksample touse

	qui count if `touse'
	if r(N) == 0 error 2000

if "`bmi'" !="" {

qui gen friskcvd30=. if `touse'
qui gen friskcvd30opt=. if `touse'
qui gen friskcvd30nor=. if `touse'

label variable friskcvd30 "General 30-year CVD risk (BMI)"
label variable friskcvd30opt "General 30-year CVD risk (optimal, BMI)"
label variable friskcvd30nor "General 30-year CVD risk (normal, BMI)"

qui gen hriskcvd30=. if `touse'
qui gen hriskcvd30opt=. if `touse'
qui gen hriskcvd30nor=. if `touse'

label variable hriskcvd30 "Hard 30-year CVD risk (BMI)"
label variable hriskcvd30opt "Hard 30-year CVD risk (optimal, BMI)"
label variable hriskcvd30nor "Hard 30-year CVD risk (normal, BMI)"

// execute program computing estimator and its covariance matrix in Mata

markout `touse'

mata: RISKF("`male'","`age'","`sbp'","`smoke'","`trtbp'","`bmi'","`diab'","`touse'")
mata: RISKH("`male'","`age'","`sbp'","`smoke'","`trtbp'","`bmi'","`diab'","`touse'")

qui replace friskcvd30=. if friskcvd30==0
qui replace friskcvd30opt=. if friskcvd30opt==0
qui replace friskcvd30nor=. if friskcvd30nor==0

qui replace hriskcvd30=. if hriskcvd30==0
qui replace hriskcvd30opt=. if hriskcvd30opt==0
qui replace hriskcvd30nor=. if hriskcvd30nor==0
}

if "`tcl'" !="" {

qui gen friskcvdl30=. if `touse'
qui gen friskcvdl30opt=. if `touse'
qui gen friskcvdl30nor=. if `touse'

label variable friskcvdl30 "General 30-year CVD risk (lipids)"
label variable friskcvdl30opt "General 30-year CVD risk (optimal, lipids)"
label variable friskcvdl30nor "General 30-year CVD risk (normal, lipids)"

qui gen hriskcvdl30=. if `touse'
qui gen hriskcvdl30opt=. if `touse'
qui gen hriskcvdl30nor=. if `touse'

label variable hriskcvdl30 "Hard 30-year CVD risk (lipids)"
label variable hriskcvdl30opt "Hard 30-year CVD risk (optimal, lipids)"
label variable hriskcvdl30nor "Hard 30-year CVD risk (normal, lipids)"

// execute program computing estimator and its covariance matrix in Mata

markout `touse'

mata: RISKLF("`male'","`age'","`sbp'","`tcl'","`hdl'","`smoke'","`trtbp'","`diab'","`touse'")
mata: RISKLH("`male'","`age'","`sbp'","`tcl'","`hdl'","`smoke'","`trtbp'","`diab'","`touse'")

qui replace friskcvdl30=. if friskcvdl30==0
qui replace friskcvdl30opt=. if friskcvdl30opt==0
qui replace friskcvdl30nor=. if friskcvdl30nor==0

qui replace hriskcvdl30=. if hriskcvdl30==0
qui replace hriskcvdl30opt=. if hriskcvdl30opt==0
qui replace hriskcvdl30nor=. if hriskcvdl30nor==0

}


end


clear mata



//------------------ Mata sub-program computing full CVD risk ------------------------------

mata

real matrix RISKF (string scalar male, string scalar age, string scalar sbp, string scalar smoke, string scalar trtbp, string scalar bmi, string scalar diab, string scalar touse)
{
real scalar N, xbeta
real vector _sex, _logage, _logsbp, _smoke, _trtbp, _logbmi, _newdiab, _inputvars,  riskvector, riskvectoropt, riskvectornor

// create matrix which is a view on the current Stata dataset

st_view(inputvars=., ., (tokens(male), tokens(age), tokens(sbp), tokens(smoke), tokens(trtbp), tokens(bmi), tokens(diab)), touse)

// inputvars[1,.]

// Transform inputed values: sex, logage, logsbp, smoke, trtbp, logbmi, newdiab

_sex=inputvars[.,1]
_logage=ln(inputvars[.,2])
_logsbp=ln(inputvars[.,3])
_smoke=inputvars[.,4]
_trtbp=inputvars[.,5]
_logbmi=ln(inputvars[.,6])
_newdiab=inputvars[.,7]

_inputvars=_sex,_logage,_logsbp,_smoke,_trtbp,_logbmi,_newdiab

//_inputvars[1,.]

// Load initial data (it's loaded in pieces due to Stata memory limitation)

// Load data for normal CVD risk

#delimit;

E1=(0.999869108564727\0.999869108564727\0.999869108564727\0.999738123610254\0.999607107747616\0.999607107747616\0.999475982600045\0.999344828842231\0.999344828842231\0.999344828842231\0.999213502103984\0.999213502103984\0.999082126266255\0.999082126266255\0.999082126266255\0.998950584114242\0.998818988324646\0.998687332389285\0.998555594029084\0.998555594029084\0.998555594029084\0.998423701998577\0.998291772647985\0.998159831430531\0.998027809366268\0.997895779523368\0.997763709100678\0.997763709100678\0.997763709100678\0.997631417994680\
0.997499094934377\0.997366695807784\0.997234283848945\0.997234283848945\0.997101829272670\0.997101829272670\0.996969309637852\0.996836724236051\0.996836724236051\0.996836724236051\0.996836724236051\0.996703641327874\0.996570530993559\0.996437396327449\0.996304176062190\0.996170890322281\0.996170890322281\0.996037413507343\0.995903656924018\0.995769870739940\0.995635999776199\0.995502093950072\0.995368060928547\0.995233925633727\0.995099585055648\0.994965129714331\0.994830594955310\0.994830594955311\0.994830594955310\0.994695815381768\
0.994560894637525\0.994425700022916\0.994290350149534\0.994290350149534\0.994290350149534\0.994154803464763\0.994019251351317\0.993883423612873\0.993747573263051\0.993611656940137\0.993475682780310\0.993339264777366\0.993202844765638\0.993066393233372\0.992929943486299\0.992793403613626\0.992656706115493\0.992519985706665\0.992383248774716\0.992246486480539\0.992109704132654\0.991972926262601\0.991836108850300\0.991699274582816\0.991562381672870\0.991425397463937\0.991288393632163\0.991151355914040\0.991014213585958\0.990876980437205\
0.990739733594364\0.990602449358768\0.990465052343700\0.990327559709196\0.990189911033584\0.990052192742558\0.989914191132021\0.989776141654459\0.989638048681040\0.989499969908606\0.989361758012742\0.989223437382275\0.989085036850462\0.988946584553709\0.988808127678336\0.988669651453813\0.988531071015477\0.988392390205740\0.988253706825485\0.988253706825485\0.988115012918886\0.987976287845952\0.987837559104631\0.987698717657703\0.987559822309340\0.987420867744568\0.987281891378805\0.987142782577555\0.987003660920324\0.986864544171682\
0.986725399225593\0.986586233873323\0.986447018411335\0.986307787496305\0.986168149234039\0.986028505497682\0.985888804367470\0.985888804367470\0.985749046022420\0.985609269468643\0.985609269468643\0.985469166714393\0.985469166714393\0.985469166714393\0.985328869737533\0.985188428360835\0.985047856846056\0.984907162120426\0.984766124566172\0.984766124566172\0.984624976524285\0.984483809131706\0.984342572640132\0.984201348815587\0.984060109796905\0.983918851797832\0.983777609796560\0.983777609796560\0.983636275411492\0.983494886744415\
0.983353350770222\0.983353350770221\0.983211652048527\0.983069943646985\0.982928179939492\0.982786236972396\0.982644182314363\0.982502080011113\0.982359867210242\0.982359867210242\0.982359867210242\0.982359867210242\0.982359867210242\0.982217219016123\0.982074529377644\0.981931785837838\0.981788781180028\0.981645763947467\0.981502739437257\0.981359518092981\0.981216271700744\0.981072927628873\0.980929578876383\0.980929578876383\0.980929578876383\0.980929578876383\0.980786058450219\0.980786058450219\0.980642493219396\0.980498838216246\
0.980355136461751\0.980211404932738\0.980067646153717\0.979923861671150\0.979779938707741\0.979635995726231\0.979635995726231\0.979491939474328\0.979347811460702\0.979203672699496\0.979059479586765\0.978915260563687\0.978771033976426\0.978626729707501\0.978482339835731\0.978337915294481\0.978193415991480\0.978048864349701\0.977904096167960\0.977904096167960\0.977904096167960\0.977759217638205\0.977614285541339\0.977614285541339\0.977469148161171\0.977469148161171\0.977323731853096\0.977323731853096\0.977178235315759\0.977032734527098\
0.977032734527098\0.976886972828746\0.976741189359697\0.976595334015378\0.976595334015378\0.976449439522738\0.976303261105114\0.976157062486451\0.976157062486451\0.976010733183493\0.975864345082117\0.975717936233427\0.975571328862367\0.975424656888729\0.975277980686276\0.975277980686277\0.975130816154825\0.974983620921294\0.974836383391717\0.974688983629767\0.974541558547654\0.974394004872218\0.974394004872218\0.974246335952206\0.974098683920345\0.973951002721556\0.973803286211087\0.973803286211087\0.973655447550031\0.973507600344844);
E2=(0.973507600344844\0.973359724131992\0.973359724131992\0.973359724131992\0.973211616901719\0.973063474801535\0.972915300575951\0.972767068191456\0.972618728997952\0.972618728997952\0.972618728997952\0.972469806113277\0.972320899575093\0.972171950371437\0.972022990985574\0.972022990985574\0.971873886251537\0.971873886251537\0.971724758042026\0.971575594239746\0.971426321382558\0.971277044588067\0.971127710814117\0.970978385569985\0.970829032740622\0.970679543273827\0.970529911145416\0.970380238776784\0.970230543083815\0.970080861129551\
0.969931177426737\0.969931177426737\0.969781282484317\0.969631322761476\0.969481250006489\0.969481250006489\0.969481250006489\0.969331030591666\0.969180689773162\0.969030299221564\0.968879912845723\0.968729360926756\0.968578776128046\0.968578776128046\0.968428031145163\0.968428031145163\0.968277140750687\0.968126134465821\0.967975013760607\0.967823816411104\0.967823816411104\0.967672262738218\0.967520690823504\0.967369105822467\0.967217480116648\0.967065851157483\0.966914120981686\0.966914120981686\0.966761863704171\0.966609460809250\
0.966457056271614\0.966304659050779\0.966304659050779\0.966152180163400\0.965999680097017\0.965847153527354\0.965694639813834\0.965542126988996\0.965389576766292\0.965236911194147\0.965084221655889\0.964931493497247\0.964778709600017\0.964625653622554\0.964472597202417\0.964319465961482\0.964319465961482\0.964319465961482\0.964166091159091\0.964012708405836\0.964012708405836\0.963859308631425\0.963859308631425\0.963705752933986\0.963552118916218\0.963398494835063\0.963244834361882\0.963244834361882\0.963090988843440\0.962937118795003\
0.962783228766694\0.962629240767920\0.962475201089487\0.962475201089487\0.962321109446693\0.962321109446693\0.962166927038474\0.962012621310525\0.962012621310525\0.961858275081212\0.961858275081212\0.961858275081212\0.961703772034570\0.961549241903299\0.961394630711029\0.961394630711029\0.961239947207787\0.961085265256698\0.960930578970880\0.960775824866686\0.960621074517781\0.960466329663195\0.960466329663195\0.960311516827595\0.960311516827595\0.960156498150627\0.960001457487159\0.959846314225438\0.959846314225438\0.959691122545839\
0.959535837959684\0.959380485400166\0.959380485400166\0.959380485400166\0.959224965924445\0.959069422069389\0.959069422069389\0.958913783173108\0.958758013225192\0.958602173670155\0.958446246639677\0.958290282980079\0.958290282980079\0.958134078798200\0.958134078798200\0.957977785962724\0.957977785962724\0.957821347386736\0.957664906265072\0.957508414322756\0.957351905821059\0.957195406188509\0.957195406188509\0.957038602251106\0.957038602251106\0.956881664605338\0.956881664605338\0.956724563741573\0.956567172716970\0.956567172716970\
0.956567172716970\0.956409521771439\0.956251874276185\0.956251874276185\0.956093962913974\0.955936034554932\0.955778020276620\0.955620010371405\0.955461859909257\0.955303543263354\0.955303543263354\0.955303543263354\0.955145030803224\0.955145030803224\0.954986319675397\0.954827567291271\0.954668816919484\0.954509999208688\0.954350797243306\0.954350797243306\0.954191531122090\0.954032232197082\0.954032232197082\0.954032232197082\0.953872676164256\0.953713058440129\0.953713058440129\0.953713058440129\0.953553313867478\0.953393297266205\
0.953393297266205\0.953393297266205\0.953233049358777\0.953072772627078\0.952912436627577\0.952752081867194\0.952591679242006\0.952431261228945\0.952270806021948\0.952270806021948\0.952110202258055\0.952110202258055\0.951949577198083\0.951788926422067\0.951628240327199\0.951628240327199\0.951467134942401\0.951467134942401\0.951467134942401\0.951305805228143\0.951144471764533\0.951144471764533\0.950983068039053\0.950821575881944\0.950821575881944\0.950660043903181\0.950498112016868\0.950336145502478\0.950174063804150\0.950011697575961\
0.949849190149732\0.949686520216158\0.949686520216158\0.949523813552735\0.949361111630044\0.949198309106864\0.949035328538795\0.948872247459750\0.948709117825486\0.948545956020237\0.948382769174259\0.948219541742828\0.948056099306325\0.947892634825951\0.947892634825951\0.947729042289161\0.947565288629925\0.947565288629925\0.947565288629925\0.947400956859161\0.947400956859161\0.947400956859161\0.947400956859161\0.947236067939443\0.947070880511413\0.946905445948185\0.946739947659505\0.946574022178144\0.946408119116216\0.946242219424853);
E3=(0.946076300909200\0.945910345645219\0.945744320472780\0.945744320472779\0.945744320472780\0.945578158836632\0.945411911174975\0.945245600683840\0.945079259137906\0.944912860907046\0.944746432114502\0.944746432114502\0.944579822846921\0.944413224061709\0.944413224061709\0.944246540013564\0.944079875238990\0.943913166965395\0.943913166965395\0.943913166965395\0.943746251594519\0.943746251594519\0.943579288883002\0.943579288883002\0.943412275848836\0.943244674495541\0.943244674495541\0.943244674495541\0.943076753516724\0.942908677185813\
0.942740400134383\0.942572129723925\0.942572129723925\0.942403623434889\0.942234992321976\0.942066357643851\0.941897696941103\0.941729044779576\0.941729044779576\0.941560068042074\0.941391070018767\0.941222083858334\0.941222083858334\0.941052984415884\0.940883882761529\0.940714699097132\0.940714699097132\0.940714699097132\0.940714699097132\0.940544499029571\0.940374216009788\0.940203716487943\0.940033196785512\0.939862682073118\0.939862682073117\0.939692000482876\0.939521306601462\0.939350571570724\0.939179808089844\0.939009047830690\
0.938838206218971\0.938838206218971\0.938667274475401\0.938496268675353\0.938496268675353\0.938325190514027\0.938154114350794\0.937982984557496\0.937982984557496\0.937982984557496\0.937811652300415\0.937640250399974\0.937468733369566\0.937468733369567\0.937297119393368\0.937125365433024\0.937125365433024\0.936953605699599\0.936953605699599\0.936781601037009\0.936609512058666\0.936609512058667\0.936437136092078\0.936437136092078\0.936437136092078\0.936437136092078\0.936264557636486\0.936091924361623\0.936091924361623\0.935919172607506\
0.935746436110172\0.935746436110173\0.935573613824872\0.935573613824872\0.935400565377945\0.935400565377945\0.935227351608029\0.935227351608028\0.935054038749162\0.934880659546370\0.934880659546370\0.934880659546370\0.934880659546370\0.934880659546370\0.934706876611052\0.934706876611052\0.934532859989287\0.934532859989287\0.934358734211608\0.934358734211608\0.934358734211608\0.934184354664466\0.934009882933250\0.933835382394768\0.933660654528032\0.933660654528032\0.933485350402887\0.933309961157497\0.933134453512749\0.932958948134548\
0.932783420606540\0.932783420606540\0.932607739278437\0.932432013400721\0.932256147368864\0.932080261137612\0.932080261137612\0.931904331363501\0.931728385788221\0.931552393846453\0.931376322619808\0.931200187146638\0.931023659407005\0.930847118805307\0.930670521887484\0.930493841401441\0.930317029372259\0.930140232825996\0.929963425420261\0.929963425420261\0.929786195970018\0.929608889953619\0.929431467344100\0.929254052251842\0.929254052251842\0.929076589498128\0.928899093696398\0.928721593783818\0.928543900693351\0.928366156427743\
0.928188390634310\0.928188390634310\0.928188390634310\0.928188390634310\0.928010131893874\0.928010131893874\0.928010131893874\0.928010131893874\0.927830588072102\0.927830588072102\0.927650892091601\0.927650892091601\0.927650892091601\0.927470961097792\0.927290983519369\0.927110954114506\0.926930927532435\0.926750846069834\0.926570597050613\0.926390265147798\0.926390265147798\0.926209770848510\0.926029164222502\0.926029164222502\0.925848399810760\0.925667631143137\0.925486886252107\0.925306155405553\0.925125446694282\0.924944650409345\
0.924763369574859\0.924582103224209\0.924582103224209\0.924400613992788\0.924219091337694\0.924037406975152\0.923855557496628\0.923855557496628\0.923673020506378\0.923490450289381\0.923307831688782\0.923125149353424\0.922942218919948\0.922759257886312\0.922759257886311\0.922575984856692\0.922392729138328\0.922209458852062\0.922026160169428\0.921842856820928\0.921659544298775\0.921476199649458\0.921292850152087\0.921292850152087\0.921109420098151\0.920925905361086\0.920925905361086\0.920742128331732\0.920558154085507\0.920374185212873\
0.920190074619921\0.920190074619922\0.920005972633612\0.919821708353122\0.919821708353122\0.919821708353122\0.919637061257534\0.919637061257534\0.919637061257535\0.919637061257534\0.919637061257534\0.919451570068011\0.919266038866714\0.919080214160695\0.919080214160695\0.918894057055487\0.918894057055487\0.918894057055487\0.918894057055487\0.918707250624897\0.918707250624897\0.918707250624897\0.918520289435702\0.918333225970073\0.918333225970073\0.918333225970073\0.918333225970073\0.918145740895012\0.918145740895012\0.917958045066648);
E4=(0.917770370249574\0.917770370249574\0.917581672082095\0.917581672082096\0.917392913530521\0.917204144128674\0.917015325545231\0.916826479553355\0.916637654238575\0.916637654238576\0.916637654238575\0.916637654238575\0.916637654238576\0.916448485993242\0.916448485993242\0.916448485993242\0.916448485993242\0.916448485993242\0.916258882602811\0.916069280165810\0.916069280165811\0.916069280165810\0.915878923996829\0.915688572617570\0.915498220569037\0.915307810708917\0.915117359136629\0.915117359136629\0.914926783105150\0.914926783105150\
0.914736147373509\0.914736147373509\0.914736147373509\0.914545132502046\0.914545132502046\0.914353921327492\0.914162668659703\0.913971390107541\0.913971390107541\0.913779925098224\0.913779925098224\0.913588358121148\0.913588358121148\0.913396735325913\0.913205105044346\0.913013421990592\0.913013421990592\0.913013421990592\0.912821613123053\0.912629782118192\0.912437929953702\0.912246072516910\0.912246072516910\0.912054101676594\0.911862077476812\0.911669969074325\0.911477832933686\0.911285666921041\0.911285666921041\0.911093354034702\
0.911093354034702\0.910900460748471\0.910707429095911\0.910707429095911\0.910514316413658\0.910514316413658\0.910514316413658\0.910320766157599\0.910127067487421\0.909932884619974\0.909738717351485\0.909544393558840\0.909544393558840\0.909544393558840\0.909544393558840\0.909349961832634\0.909155511028704\0.908960961199699\0.908960961199699\0.908766381557314\0.908571816235898\0.908571816235898\0.908377215807572\0.908377215807572\0.908377215807572\0.908377215807572\0.908182282061151\0.907987210423774\0.907792005352507\0.907792005352507\
0.907596482443454\0.907400942347424\0.907205423910710\0.907009689370304\0.906813759232372\0.906617843130188\0.906421866171257\0.906225821473271\0.906225821473271\0.906225821473271\0.906029548384792\0.906029548384792\0.905833156960828\0.905636747509218\0.905440340096624\0.905243922640983\0.905047478562650\0.904850995110113\0.904654507792354\0.904457889429585\0.904261261638310\0.904261261638309\0.904064169989894\0.904064169989894\0.903866852357655\0.903866852357655\0.903669410960330\0.903669410960330\0.903471641278071\0.903273855511004\
0.903076000836677\0.903076000836677\0.903076000836678\0.902878041861201\0.902680016695893\0.902481943187683\0.902283838643601\0.902283838643601\0.902085671453439\0.901887395645964\0.901689094494440\0.901490827155985\0.901292542642492\0.901094189229688\0.901094189229688\0.901094189229688\0.901094189229688\0.901094189229688\0.900894676811393\0.900695156678769\0.900695156678769\0.900695156678769\0.900494581355491\0.900494581355491\0.900293988657176\0.900093428743137\0.899892850977497\0.899892850977498\0.899892850977497\0.899692054229204\
0.899490782574674\0.899289438673322\0.899088073886084\0.898886696959037\0.898685184035910\0.898483582450614\0.898281866656168\0.898080176550853\0.898080176550853\0.897878433916939\0.897878433916939\0.897676080393873\0.897676080393873\0.897473657799038\0.897271174158623\0.897271174158623\0.897271174158623\0.897271174158622\0.897068478162465\0.896865789207261\0.896865789207261\0.896662944363780\0.896459997254760\0.896256970280460\0.896053908649170\0.895850821584265\0.895850821584265\0.895850821584265\0.895647456030898\0.895443972631846\
0.895240363276796\0.895036655456051\0.894832948782125\0.894832948782125\0.894629105037506\0.894425268427283\0.894425268427283\0.894425268427283\0.894221143708050\0.894221143708050\0.894016768700476\0.893812347931776\0.893607855240777\0.893403374844241\0.893403374844241\0.893198787885899\0.893198787885899\0.892994195564994\0.892789520919526\0.892584704589431\0.892379736071846\0.892379736071846\0.892174516874986\0.891969151671544\0.891763520130245\0.891763520130245\0.891557588781600\0.891351623484210\0.891145606007355\0.890939548800815\
0.890733374168164\0.890527195972940\0.890527195972941\0.890527195972940\0.890527195972940\0.890320488990133\0.890320488990133\0.890113522464911\0.889906332604333\0.889699130297180\0.889491859654277\0.889284526739487\0.889077170088466\0.889077170088466\0.888869324264266\0.888869324264266\0.888661429530998\0.888453497618032\0.888453497618032\0.888245417393924\0.888037247260640\0.888037247260640\0.887828898972975\0.887620421173760\0.887620421173760\0.887411662368718\0.887411662368717\0.887202811946717\0.887202811946717\0.886993868433346);
E5=(0.886993868433346\0.886784663473047\0.886784663473047\0.886575237289968\0.886365740610301\0.886365740610301\0.886156165902640\0.885946461518301\0.885946461518301\0.885736506666694\0.885526558929437\0.885316597350439\0.885106658265826\0.884896680182418\0.884686574542743\0.884686574542743\0.884475991386741\0.884265118660875\0.884054144597343\0.884054144597343\0.883842957672929\0.883631512162798\0.883419624948073\0.883207684153032\0.882995724721651\0.882783760597658\0.882783760597658\0.882783760597658\0.882783760597658\0.882783760597658\
0.882783760597658\0.882571097037874\0.882358417265080\0.882358417265080\0.882145709641802\0.881932911491781\0.881719974095919\0.881506955934083\0.881506955934083\0.881506955934083\0.881293634303050\0.881293634303050\0.881080217985702\0.880866811236585\0.880653415454246\0.880653415454246\0.880653415454247\0.880439621953710\0.880225747817598\0.880011797634808\0.880011797634808\0.880011797634808\0.879797647428018\0.879797647428018\0.879583449039404\0.879369274251383\0.879155111157026\0.879155111157026\0.878940780633317\0.878940780633317\
0.878940780633317\0.878940780633317\0.878726304074691\0.878726304074691\0.878511588356346\0.878511588356346\0.878296641322989\0.878081495113038\0.877866352427736\0.877651251115045\0.877651251115045\0.877436016662925\0.877220743323482\0.877005315124943\0.876789884995100\0.876574441526077\0.876359018557289\0.876359018557288\0.876143516588475\0.875927874848740\0.875927874848740\0.875712007663958\0.875495829726444\0.875279603602283\0.875279603602283\0.875063256601253\0.875063256601253\0.874846855980563\0.874630465531000\0.874414022656760\
0.874197611923066\0.874197611923066\0.873980831341721\0.873764100061325\0.873764100061325\0.873547222974077\0.873547222974078\0.873547222974078\0.873329626292952\0.873329626292952\0.873329626292952\0.873111885494356\0.872894063124732\0.872676134353552\0.872458098006720\0.872239906901864\0.872021646459213\0.872021646459213\0.871803134255671\0.871803134255671\0.871803134255671\0.871584327000165\0.871584327000165\0.871584327000165\0.871364909119256\0.871145441615907\0.870925941252386\0.870706359631889\0.870706359631889\0.870486585188871\
0.870486585188871\0.870266683172814\0.870046718317388\0.869826657342776\0.869826657342776\0.869826657342776\0.869826657342776\0.869605772847298\0.869384925793766\0.869384925793766\0.869384925793766\0.869163880104542\0.868942879646934\0.868942879646934\0.868721796185747\0.868500545594694\0.868279040453837\0.868279040453837\0.868057182370107\0.868057182370107\0.867835267361812\0.867835267361811\0.867612982602805\0.867390629486898\0.867168186324258\0.867168186324258\0.866945328608073\0.866722502924560\0.866722502924559\0.866499662574462\
0.866276561611604\0.866053332717352\0.865830068846433\0.865606267845330\0.865382441994783\0.865158579063614\0.864934627856751\0.864710722084797\0.864486770081876\0.864486770081876\0.864486770081876\0.864262066466850\0.864036909767367\0.864036909767366\0.863811653519634\0.863586351321963\0.863361029606081\0.863361029606081\0.863135514946389\0.863135514946389\0.863135514946389\0.862909458988727\0.862909458988727\0.862909458988727\0.862683221717449\0.862683221717449\0.862683221717449\0.862456119475797\0.862228985932829\0.862001600135367\
0.861774228034333\0.861774228034333\0.861546487963908\0.861546487963908\0.861546487963908\0.861318534870280\0.861090538481091\0.860862527291299\0.860634393829357\0.860406288152326\0.860406288152326\0.860406288152326\0.860177924448274\0.859949498581782\0.859721000554280\0.859721000554280\0.859491908560172\0.859262838176271\0.859033743756790\0.858804654060891\0.858575368955137\0.858345848052403\0.858345848052403\0.858345848052403\0.858345848052403\0.858345848052403\0.858345848052403\0.858345848052403\0.858114729392155\0.857883636669681\
0.857652570394043\0.857421483310152\0.857421483310152\0.857190199613556\0.857190199613557\0.856958913904039\0.856727537019876\0.856495957119963\0.856264265650567\0.856264265650567\0.856264265650567\0.856032115826013\0.855800002719949\0.855567892708408\0.855335747454557\0.855335747454557\0.855335747454557\0.855335747454557\0.855335747454558\0.855102598493535\0.855102598493535\0.855102598493535\0.855102598493535\0.854869060681614\0.854869060681614\0.854635309392208\0.854635309392207\0.854401219854382\0.854401219854383\0.854401219854382);
E6=(0.854166660052329\0.853931974057729\0.853931974057729\0.853697010671110\0.853697010671110\0.853697010671110\0.853461803083981\0.853461803083981\0.853461803083981\0.853226123206046\0.853226123206046\0.852990147551823\0.852754032268753\0.852517762583941\0.852281466732581\0.852281466732581\0.852045037851344\0.852045037851344\0.851808580318012\0.851571771452307\0.851334657162115\0.851097460883225\0.850860251637694\0.850622895313322\0.850622895313322\0.850385373791552\0.850147814541005\0.849910303188567\0.849672442730275\0.849434301202152\
0.849434301202152\0.849434301202152\0.849195753806659\0.849195753806659\0.848957170554385\0.848957170554385\0.848718460431687\0.848479334236427\0.848479334236427\0.848479334236427\0.848239975703147\0.848000183667034\0.847760220925750\0.847520006179913\0.847279746149370\0.847279746149370\0.847279746149370\0.847279746149370\0.847279746149370\0.847038119071354\0.847038119071353\0.847038119071353\0.846795139661503\0.846552151591903\0.846552151591903\0.846308937199885\0.846065711485289\0.845822092598133\0.845822092598133\0.845822092598134\
0.845577543324219\0.845333009306716\0.845088014611671\0.845088014611671\0.844842625127957\0.844842625127957\0.844842625127957\0.844596868907644\0.844351104049668\0.844351104049668\0.844351104049668\0.844105121648641\0.843858955514615\0.843858955514615\0.843612558647887\0.843365895591231\0.843365895591231\0.843118987086256\0.842872051952526\0.842624655273729\0.842377306367796\0.842377306367796\0.842129424582218\0.841881465320426\0.841633349733917\0.841385269837213\0.841137187215830\0.841137187215830\0.840888886458160\0.840888886458159\
0.840640443017062\0.840391849074818\0.840143260876721\0.839894647482556\0.839645446709836\0.839396241497589\0.839396241497589\0.839396241497589\0.839146506821057\0.839146506821057\0.839146506821057\0.838896007640840\0.838896007640840\0.838896007640840\0.838644303902772\0.838644303902772\0.838391701338560\0.838139045194927\0.837885653954454\0.837631596474656\0.837377257998536\0.837122597167312\0.836867669365043\0.836867669365043\0.836867669365043\0.836867669365043\0.836610953584868\0.836610953584868\0.836353738806831\0.836353738806831\
0.836096243747031\0.835838788867230\0.835838788867230\0.835838788867230\0.835838788867230\0.835579453112403\0.835319684884801\0.835319684884801\0.835319684884801\0.835319684884801\0.835058371176515\0.835058371176516\0.834796360524142\0.834796360524141\0.834533355986523\0.834533355986523\0.834269696015794\0.834269696015794\0.834269696015794\0.834269696015794);

J1=(1.000000000000000\0.999874999925279\0.999749974759799\0.999749974759799\0.999749974759799\0.999624814063019\0.999624814063019\0.999624814063019\0.999499490200862\0.999374072344002\0.999374072344002\0.999248569418730\0.999248569418730\0.999123023243039\0.998997386740371\0.998997386740371\0.998997386740371\0.998997386740371\0.998997386740371\0.998871256730167\0.998745136123117\0.998745136123117\0.998745136123117\0.998745136123117\0.998745136123117\0.998745136123117\0.998745136123117\0.998618666245649\0.998492024717242\0.998492024717242\
0.998492024717242\0.998492024717242\0.998492024717242\0.998365144913641\0.998365144913641\0.998238163160071\0.998238163160071\0.998238163160071\0.998110958718562\0.997983509789249\0.997856068785719\0.997856068785719\0.997856068785719\0.997856068785719\0.997856068785719\0.997856068785719\0.997728162803535\0.997728162803535\0.997728162803535\0.997728162803535\0.997728162803535\0.997728162803535\0.997728162803535\0.997728162803535\0.997728162803535\0.997728162803535\0.997728162803535\0.997598893891575\0.997469610578858\0.997469610578858\
0.997469610578858\0.997469610578858\0.997469610578858\0.997339669327829\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\
0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997209721540522\0.997076012382359\0.997076012382359\0.997076012382359\0.997076012382359\0.997076012382359\0.997076012382359\0.997076012382359\0.997076012382359\0.997076012382359\0.997076012382359\0.997076012382359\
0.997076012382359\0.997076012382359\0.997076012382359\0.997076012382359\0.997076012382359\0.997076012382359\0.997076012382359\0.996941193341222\0.996941193341222\0.996941193341222\0.996806244147050\0.996806244147050\0.996670918946664\0.996535466441059\0.996535466441059\0.996535466441059\0.996535466441059\0.996535466441059\0.996535466441059\0.996399076293597\0.996399076293597\0.996399076293597\0.996399076293597\0.996399076293597\0.996399076293597\0.996399076293597\0.996399076293597\0.996262427661643\0.996262427661643\0.996262427661643\
0.996262427661643\0.996125496615110\0.996125496615110\0.996125496615110\0.996125496615110\0.996125496615110\0.996125496615110\0.996125496615110\0.996125496615110\0.995987653348361\0.995849808215122\0.995711894055637\0.995573855258422\0.995573855258422\0.995573855258422\0.995573855258422\0.995573855258422\0.995573855258422\0.995573855258422\0.995573855258422\0.995573855258422\0.995573855258422\0.995573855258422\0.995434919915833\0.995295967956718\0.995156976694060\0.995156976694060\0.995017843135763\0.995017843135763\0.995017843135763\
0.995017843135763\0.995017843135763\0.995017843135763\0.995017843135763\0.995017843135763\0.995017843135763\0.994878162415113\0.994878162415113\0.994878162415113\0.994878162415113\0.994878162415113\0.994878162415113\0.994878162415113\0.994878162415113\0.994878162415113\0.994878162415113\0.994878162415113\0.994878162415113\0.994878162415113\0.994737447278356\0.994596746095655\0.994596746095655\0.994596746095655\0.994455713487270\0.994455713487270\0.994314623000990\0.994314623000990\0.994173217376508\0.994173217376508\0.994173217376508\
0.994031466120154\0.994031466120154\0.994031466120154\0.994031466120154\0.993889472916219\0.993889472916219\0.993889472916219\0.993889472916219\0.993747047555321\0.993747047555321\0.993747047555321\0.993747047555321\0.993747047555321\0.993747047555321\0.993747047555321\0.993604074592034\0.993604074592034\0.993604074592034\0.993604074592034\0.993604074592034\0.993604074592034\0.993604074592034\0.993460131547987\0.993460131547987\0.993460131547987\0.993460131547987\0.993460131547987\0.993316068706432\0.993316068706432\0.993316068706432);
J2=(0.993171893774472\0.993171893774472\0.993027652598441\0.992883222236069\0.992883222236069\0.992883222236069\0.992883222236069\0.992883222236069\0.992883222236069\0.992738306041608\0.992592880224925\0.992592880224925\0.992592880224925\0.992592880224925\0.992592880224925\0.992447191704021\0.992447191704021\0.992301390716518\0.992301390716518\0.992301390716518\0.992301390716518\0.992301390716518\0.992301390716518\0.992301390716518\0.992301390716518\0.992301390716518\0.992301390716518\0.992301390716518\0.992301390716518\0.992301390716518\
0.992301390716518\0.992154765161896\0.992154765161896\0.992154765161896\0.992154765161896\0.992007901712901\0.991860984449752\0.991860984449752\0.991860984449752\0.991860984449752\0.991860984449752\0.991860984449752\0.991860984449752\0.991713502037289\0.991713502037289\0.991565821218191\0.991565821218191\0.991565821218191\0.991565821218191\0.991565821218191\0.991417685304608\0.991417685304608\0.991417685304608\0.991417685304608\0.991417685304608\0.991417685304608\0.991417685304608\0.991268920479898\0.991268920479898\0.991268920479898\
0.991268920479898\0.991268920479898\0.991119606551472\0.991119606551472\0.991119606551472\0.991119606551472\0.991119606551472\0.991119606551472\0.991119606551472\0.991119606551472\0.991119606551472\0.991119606551472\0.991119606551472\0.991119606551472\0.991119606551472\0.991119606551472\0.990969335120654\0.990818995577878\0.990818995577878\0.990818995577878\0.990668469736001\0.990668469736001\0.990517911061998\0.990517911061998\0.990517911061998\0.990517911061998\0.990517911061998\0.990366962714409\0.990366962714409\0.990366962714409\
0.990366962714409\0.990366962714409\0.990366962714409\0.990215666076276\0.990215666076276\0.990064267775677\0.990064267775677\0.990064267775677\0.989912695203268\0.989912695203268\0.989761073478035\0.989609461419450\0.989609461419450\0.989609461419450\0.989609461419450\0.989457535334001\0.989457535334001\0.989457535334001\0.989457535334001\0.989457535334001\0.989457535334001\0.989457535334001\0.989305319775663\0.989305319775663\0.989153024957103\0.989153024957103\0.989153024957103\0.989153024957103\0.989000470974402\0.989000470974402\
0.989000470974402\0.989000470974402\0.988847593883316\0.988694612560821\0.988694612560821\0.988694612560821\0.988541472767695\0.988541472767695\0.988541472767695\0.988541472767695\0.988541472767695\0.988541472767695\0.988387775436182\0.988387775436182\0.988233876011081\0.988233876011081\0.988079805661552\0.988079805661552\0.988079805661552\0.988079805661552\0.988079805661552\0.988079805661552\0.987925306390708\0.987925306390708\0.987770672122576\0.987770672122576\0.987615840452036\0.987615840452036\0.987615840452036\0.987460338805120\
0.987304769234185\0.987304769234185\0.987304769234185\0.987148887094422\0.987148887094422\0.987148887094422\0.987148887094422\0.987148887094422\0.987148887094422\0.987148887094422\0.986992409002289\0.986835925995771\0.986835925995771\0.986679214388418\0.986679214388418\0.986679214388418\0.986679214388418\0.986679214388418\0.986679214388418\0.986521493064814\0.986521493064814\0.986521493064814\0.986363479965186\0.986205372304753\0.986205372304753\0.986205372304753\0.986047066493411\0.985888677702618\0.985888677702618\0.985888677702618\
0.985730091200487\0.985571407266251\0.985571407266251\0.985571407266251\0.985571407266251\0.985571407266251\0.985571407266251\0.985571407266251\0.985571407266251\0.985412287358133\0.985412287358133\0.985253038362047\0.985253038362047\0.985253038362047\0.985253038362047\0.985093417938402\0.985093417938402\0.984933574358521\0.984773697378992\0.984773697378992\0.984773697378992\0.984613736982138\0.984613736982138\0.984613736982138\0.984453547670246\0.984453547670246\0.984453547670246\0.984453547670246\0.984453547670246\0.984453547670246\
0.984453547670246\0.984453547670246\0.984291646102087\0.984291646102087\0.984291646102087\0.984291646102087\0.984291646102087\0.984291646102087\0.984291646102087\0.984291646102087\0.984291646102087\0.984291646102087\0.984291646102087\0.984291646102087\0.984128608959539\0.984128608959539\0.984128608959539\0.983965190207230\0.983801611852957\0.983801611852957\0.983637749259413\0.983473825709107\0.983309856792326\0.983309856792326\0.983309856792326\0.983309856792326\0.983309856792326\0.983309856792326\0.983309856792326\0.983309856792326);
J3=(0.983309856792326\0.983309856792326\0.983309856792326\0.983144714435010\0.982979561540553\0.982979561540553\0.982979561540553\0.982979561540553\0.982979561540553\0.982979561540553\0.982979561540553\0.982813678949649\0.982813678949649\0.982813678949649\0.982647674058965\0.982647674058965\0.982647674058965\0.982647674058965\0.982481451283708\0.982315219066056\0.982315219066056\0.982148883311126\0.982148883311126\0.981982500152994\0.981982500152994\0.981982500152994\0.981815744741301\0.981648749572087\0.981648749572087\0.981648749572087\
0.981648749572087\0.981648749572087\0.981480976760355\0.981480976760355\0.981480976760355\0.981480976760355\0.981480976760355\0.981480976760355\0.981312697811611\0.981312697811611\0.981312697811611\0.981312697811611\0.981144025358919\0.981144025358919\0.981144025358919\0.981144025358919\0.980974895346117\0.980805654738607\0.980636378506403\0.980636378506403\0.980636378506403\0.980636378506403\0.980636378506403\0.980636378506403\0.980466421656547\0.980466421656547\0.980466421656547\0.980466421656547\0.980466421656547\0.980466421656547\
0.980466421656547\0.980296001823715\0.980296001823714\0.980296001823715\0.980125402501540\0.980125402501540\0.980125402501539\0.980125402501539\0.979954543728665\0.979783689215987\0.979783689215987\0.979783689215987\0.979783689215987\0.979612514269731\0.979612514269731\0.979612514269731\0.979441192442454\0.979441192442454\0.979269712541791\0.979269712541791\0.979269712541791\0.979097926735592\0.979097926735592\0.978925803088803\0.978753692706232\0.978581589208028\0.978581589208028\0.978581589208028\0.978409307821821\0.978409307821821\
0.978409307821821\0.978236994621852\0.978236994621852\0.978064439054746\0.978064439054746\0.977891668614707\0.977891668614707\0.977718832831767\0.977718832831767\0.977718832831767\0.977545703465066\0.977372529529600\0.977199361668416\0.977026105845750\0.977026105845750\0.976852415180112\0.976852415180112\0.976678586742615\0.976678586742615\0.976504540092880\0.976330484677873\0.976330484677873\0.976330484677873\0.976330484677873\0.976330484677873\0.976155656215238\0.976155656215238\0.976155656215238\0.976155656215238\0.976155656215238\
0.976155656215238\0.975980387057773\0.975980387057773\0.975980387057773\0.975980387057773\0.975980387057773\0.975804750164314\0.975804750164314\0.975804750164314\0.975804750164314\0.975804750164314\0.975804750164314\0.975804750164314\0.975804750164314\0.975804750164314\0.975804750164314\0.975804750164314\0.975804750164314\0.975804750164314\0.975627656289966\0.975627656289966\0.975627656289966\0.975627656289966\0.975627656289966\0.975450148330391\0.975450148330391\0.975450148330391\0.975450148330391\0.975450148330391\0.975450148330391\
0.975450148330391\0.975272007255045\0.975093621833901\0.974915230008115\0.974915230008115\0.974735950640477\0.974556502179288\0.974376722361019\0.974376722361019\0.974196693860591\0.974196693860591\0.974016567282371\0.973836435883878\0.973836435883878\0.973836435883878\0.973836435883878\0.973836435883878\0.973836435883878\0.973836435883878\0.973836435883878\0.973655210804627\0.973655210804627\0.973655210804627\0.973473809127167\0.973473809127167\0.973473809127167\0.973473809127167\0.973473809127167\0.973473809127167\0.973473809127167\
0.973473809127167\0.973473809127167\0.973291342087169\0.973291342087169\0.973291342087169\0.973291342087169\0.973291342087169\0.973107918476482\0.973107918476482\0.973107918476482\0.973107918476482\0.973107918476482\0.973107918476482\0.973107918476482\0.972923354111095\0.972923354111095\0.972923354111095\0.972923354111095\0.972923354111095\0.972923354111095\0.972923354111095\0.972923354111095\0.972923354111095\0.972738277669984\0.972738277669984\0.972738277669984\0.972552568896175\0.972552568896175\0.972552568896175\0.972552568896175\
0.972552568896175\0.972366528727863\0.972366528727863\0.972366528727863\0.972180222076809\0.971993777542034\0.971993777542034\0.971807034858290\0.971620292087762\0.971433371226307\0.971246275245988\0.971246275245988\0.971246275245988\0.971246275245988\0.971058414703048\0.971058414703048\0.970870241755847\0.970681764441126\0.970493191106459\0.970493191106459\0.970304368417810\0.970115428189880\0.970115428189880\0.970115428189880\0.969926309172176\0.969737177004050\0.969547971640692\0.969547971640692\0.969358372606403\0.969358372606403);
J4=(0.969358372606403\0.969168410085766\0.969168410085766\0.968977769880800\0.968977769880800\0.968977769880800\0.968977769880800\0.968977769880800\0.968977769880800\0.968786741466653\0.968595736363170\0.968404695049884\0.968213645019446\0.968213645019446\0.968022524412690\0.967831404438472\0.967640143296484\0.967448873013752\0.967448873013752\0.967448873013752\0.967257216949463\0.967065402608098\0.967065402608098\0.967065402608098\0.967065402608098\0.967065402608098\0.967065402608098\0.966873153816590\0.966873153816590\0.966680707114197\
0.966680707114197\0.966488246069497\0.966295791796072\0.966295791796072\0.966103036407682\0.966103036407681\0.966103036407682\0.966103036407681\0.965909862046296\0.965909862046296\0.965716459104418\0.965716459104418\0.965523000310255\0.965523000310255\0.965523000310255\0.965523000310255\0.965329396559593\0.965135795291130\0.965135795291130\0.965135795291130\0.965135795291130\0.965135795291130\0.964941835466798\0.964941835466798\0.964941835466798\0.964941835466798\0.964941835466798\0.964941835466798\0.964747447814169\0.964747447814169\
0.964552916404511\0.964552916404511\0.964552916404511\0.964357637305898\0.964357637305898\0.964162177705644\0.963966441715593\0.963966441715593\0.963966441715593\0.963966441715593\0.963966441715593\0.963966441715593\0.963769584559161\0.963572713453857\0.963375872794764\0.963375872794764\0.963375872794764\0.963375872794764\0.963178767586686\0.963178767586686\0.963178767586686\0.962981636944804\0.962981636944804\0.962784321772847\0.962587003595379\0.962389705849871\0.962389705849871\0.962389705849871\0.962389705849871\0.962191679287199\
0.962191679287200\0.962191679287199\0.962191679287199\0.962191679287200\0.962191679287199\0.962191679287200\0.962191679287199\0.962191679287199\0.961992445715018\0.961793151776392\0.961793151776392\0.961593797354461\0.961593797354461\0.961593797354461\0.961593797354461\0.961593797354461\0.961593797354461\0.961593797354461\0.961593797354461\0.961593797354461\0.961593797354461\0.961393616702066\0.961393616702066\0.961192949519621\0.961192949519621\0.960992191555151\0.960992191555151\0.960791334504964\0.960791334504964\0.960791334504964\
0.960791334504964\0.960590127515377\0.960388894330778\0.960388894330778\0.960388894330778\0.960388894330778\0.960388894330778\0.960187236155301\0.960187236155301\0.960187236155301\0.960187236155301\0.960187236155301\0.960187236155301\0.960187236155301\0.959985072543118\0.959782569545601\0.959579962568045\0.959377279554651\0.959377279554651\0.959377279554651\0.959174191954156\0.958970732306294\0.958970732306294\0.958767228469952\0.958767228469952\0.958767228469952\0.958767228469952\0.958563579294544\0.958359953463964\0.958359953463964\
0.958359953463964\0.958359953463964\0.958359953463964\0.958359953463964\0.958359953463964\0.958359953463964\0.958359953463964\0.958359953463964\0.958155188213799\0.958155188213799\0.957949897780006\0.957949897780006\0.957744604639129\0.957744604639129\0.957744604639129\0.957538903063436\0.957333175715266\0.957127483568540\0.957127483568541\0.957127483568541\0.956921676818985\0.956921676818985\0.956921676818985\0.956921676818985\0.956921676818985\0.956921676818985\0.956715017844550\0.956508326963985\0.956508326963985\0.956508326963985\
0.956508326963985\0.956508326963985\0.956508326963985\0.956301061188364\0.956301061188364\0.956301061188364\0.956093625358714\0.955885925288270\0.955885925288270\0.955677982052025\0.955677982052025\0.955677982052025\0.955677982052025\0.955677982052025\0.955469603395102\0.955469603395102\0.955261167652013\0.955261167652013\0.955261167652013\0.955261167652013\0.955261167652013\0.955052242611016\0.955052242611016\0.955052242611016\0.955052242611016\0.954842663173657\0.954842663173657\0.954842663173657\0.954842663173657\0.954842663173657\
0.954842663173657\0.954842663173657\0.954632487044142\0.954422063726574\0.954211625451043\0.954211625451043\0.954000999205141\0.954000999205141\0.954000999205141\0.954000999205141\0.954000999205141\0.954000999205141\0.954000999205141\0.953789320899351\0.953789320899351\0.953577534503218\0.953577534503218\0.953577534503218\0.953365598747052\0.953365598747052\0.953365598747052\0.953153361590458\0.953153361590458\0.953153361590458\0.952940789982197\0.952940789982197\0.952727993052955\0.952727993052955\0.952515046599274\0.952515046599274);
J5=(0.952301929149243\0.952301929149243\0.952088318264917\0.952088318264917\0.952088318264917\0.951874353667158\0.951874353667158\0.951874353667158\0.951660125011556\0.951660125011556\0.951660125011556\0.951660125011555\0.951660125011556\0.951660125011556\0.951660125011556\0.951445215480848\0.951445215480848\0.951445215480848\0.951445215480848\0.951229428453214\0.951229428453214\0.951229428453214\0.951229428453214\0.951229428453214\0.951229428453214\0.951229428453214\0.951012679964095\0.950795935937874\0.950579033803043\0.950362113253201\
0.950145140081369\0.950145140081369\0.950145140081369\0.949928086592563\0.949928086592563\0.949928086592563\0.949928086592563\0.949928086592563\0.949710101096608\0.949491975074546\0.949491975074545\0.949273858516207\0.949273858516207\0.949273858516207\0.949273858516207\0.949055439619586\0.948837029043640\0.948837029043640\0.948837029043640\0.948837029043640\0.948618250373515\0.948399504946574\0.948399504946574\0.948180569563260\0.948180569563260\0.948180569563260\0.948180569563260\0.947961442488587\0.947961442488587\0.947742123769990\
0.947522736407107\0.947303362608427\0.947303362608427\0.947083957689065\0.947083957689065\0.946864168616785\0.946864168616785\0.946864168616785\0.946864168616785\0.946864168616785\0.946643928927254\0.946643928927254\0.946643928927254\0.946643928927254\0.946643928927254\0.946643928927254\0.946643928927254\0.946423066157614\0.946423066157614\0.946423066157614\0.946201846168126\0.946201846168126\0.946201846168126\0.946201846168126\0.945980086482764\0.945980086482763\0.945758098778893\0.945758098778893\0.945758098778893\0.945758098778893\
0.945758098778893\0.945535432654063\0.945535432654063\0.945535432654063\0.945312677235573\0.945312677235573\0.945089264001831\0.944865715539511\0.944865715539512\0.944641994227255\0.944418189406701\0.944418189406701\0.944418189406701\0.944418189406701\0.944418189406701\0.944418189406701\0.944418189406701\0.944193518469619\0.944193518469619\0.943968671285045\0.943743634727189\0.943743634727189\0.943518391830524\0.943292765297592\0.943292765297592\0.943292765297592\0.943292765297592\0.943292765297592\0.943066362358650\0.943066362358650\
0.942839755406199\0.942839755406200\0.942839755406199\0.942839755406200\0.942612712482589\0.942385146250573\0.942157365963273\0.942157365963273\0.942157365963273\0.941929414673893\0.941701485049577\0.941701485049577\0.941701485049577\0.941473359894618\0.941473359894618\0.941473359894618\0.941473359894618\0.941244445110600\0.941244445110600\0.941015145545997\0.941015145545997\0.940785803015591\0.940785803015591\0.940785803015591\0.940785803015591\0.940555549101437\0.940555549101437\0.940555549101437\0.940325057094620\0.940325057094620\
0.940325057094620\0.940325057094620\0.940325057094620\0.940325057094619\0.940325057094620\0.940325057094620\0.940325057094619\0.940325057094620\0.940325057094620\0.940092868587684\0.939860383804563\0.939860383804563\0.939860383804563\0.939626720768001\0.939626720768001\0.939626720768001\0.939626720768001\0.939392719404964\0.939392719404964\0.939158667854152\0.938924583816651\0.938924583816651\0.938689993195969\0.938455419725454\0.938455419725454\0.938220246741049\0.937984507304453\0.937984507304453\0.937984507304453\0.937984507304453\
0.937984507304453\0.937748030475488\0.937748030475488\0.937511521236963\0.937275028850107\0.937275028850107\0.937275028850107\0.937275028850107\0.937275028850107\0.937275028850107\0.937037926146014\0.936800845978481\0.936800845978481\0.936800845978481\0.936800845978481\0.936563266854809\0.936563266854809\0.936563266854809\0.936563266854809\0.936563266854809\0.936563266854809\0.936563266854809\0.936324802303127\0.936086357823265\0.935847741546383\0.935608439982274\0.935369151151008\0.935129741980994\0.935129741980994\0.935129741980994\
0.935129741980994\0.935129741980994\0.934890052060519\0.934890052060519\0.934650291028054\0.934650291028054\0.934650291028054\0.934650291028054\0.934650291028054\0.934409758777445\0.934168829727035\0.934168829727035\0.934168829727035\0.934168829727035\0.934168829727035\0.933927638968959\0.933686454940939\0.933445129297581\0.933203834001824\0.933203834001824\0.932961893052735\0.932719908067458\0.932477803073873\0.932477803073872\0.932235534643321\0.932235534643321\0.931992961351262\0.931992961351262\0.931750231476786\0.931507425328436);
J6=(0.931507425328436\0.931507425328436\0.931264198313636\0.931264198313636\0.931020822732180\0.930777154005075\0.930777154005075\0.930533205827668\0.930289242532862\0.930289242532862\0.930044922247188\0.930044922247188\0.930044922247188\0.930044922247188\0.930044922247188\0.929799595461934\0.929799595461934\0.929554187972832\0.929554187972832\0.929554187972832\0.929554187972832\0.929554187972832\0.929554187972832\0.929554187972832\0.929307123884548\0.929307123884548\0.929307123884548\0.929307123884548\0.929307123884548\0.929307123884548\
0.929059300491410\0.928811343338812\0.928811343338812\0.928563277074065\0.928563277074065\0.928315081097672\0.928315081097672\0.928315081097672\0.928066104556865\0.927816836487737\0.927816836487737\0.927816836487737\0.927816836487737\0.927816836487737\0.927816836487737\0.927565685278048\0.927314175585628\0.927062623407062\0.926810973312167\0.926810973312167\0.926558508301640\0.926305407396630\0.926305407396630\0.926305407396630\0.926051933674402\0.926051933674402\0.926051933674402\0.926051933674402\0.925797382959281\0.925542633759220\
0.925542633759220\0.925542633759220\0.925542633759220\0.925286993649247\0.925286993649247\0.925030852334887\0.924774682585534\0.924774682585534\0.924774682585534\0.924518251749223\0.924261831882423\0.924261831882423\0.924261831882423\0.924005098419332\0.924005098419332\0.924005098419332\0.923747652309302\0.923747652309302\0.923747652309303\0.923747652309303\0.923747652309303\0.923489860785943\0.923489860785943\0.923489860785943\0.923489860785944\0.923489860785944\0.923489860785944\0.923231108159415\0.923231108159414\0.922972260189867\
0.922972260189867\0.922972260189867\0.922972260189867\0.922972260189867\0.922972260189867\0.922972260189867\0.922712142445517\0.922451690914163\0.922451690914163\0.922190916912602\0.921929649818324\0.921929649818324\0.921667895180989\0.921405712585474\0.921405712585474\0.921142921049009\0.921142921049009\0.921142921049009\0.921142921049008\0.921142921049008\0.921142921049008\0.921142921049009\0.921142921049009\0.920875820667197\0.920608343312234\0.920340618943881\0.920340618943881\0.920071807297161\0.920071807297161\0.919802826872968\
0.919802826872969\0.919802826872968\0.919533418805403\0.919263755287570\0.918993286009764\0.918993286009764\0.918993286009764\0.918721812039330\0.918448956914244\0.918176027966571\0.918176027966571\0.917902463480191\0.917902463480191\0.917628024889568\0.917628024889567\0.917353131632463\0.917353131632463\0.917077454581598\0.916801734519996\0.916525787829497);

#delimit cr

// assemble pieces of initial data into single vectors

E=E1\E2\E3\E4\E5\E6
J=J1\J2\J3\J4\J5\J6

// Create matrix G (used below) - log{i-1}-log{i}, subset |top_left_row,top_left_column\bottom_right_row,bottom_right_column|
E1=E[|1,.\rows(E)-1,.|]
E2=E[|2,.\rows(E),.|]
G=ln(E1)-ln(E2)

// input initial values of parameters

coeffx=(0.54089\2.7684\1.69824\0.68916\0.53163\1.01325\0.79624)
coeffdx=(0.47835\3.43911\1.57043\0.98257\0.18871\-0.52377\0.51411)
mxbeta=21.8928137813482
mdxbeta=18.8268423627735

// Define new vector containing risk scores
riskvector=J(rows(_inputvars),1,.)
riskvectoropt=J(rows(_inputvars),1,.)
riskvectornor=J(rows(_inputvars),1,.)

// loop over observations and compute risk score for each one

for (i=1; i<=rows(_inputvars); i++) {

	// select observation i of input data
	Xrow=_inputvars[| i,1\i,.|]

	//--------------------- Risk score ---------------------

	// Q - make input vector into column vector
	X=Xrow'

	// R - X*beta
	xbeta=X'*coeffx
	dxbeta=X'*coeffdx

	// $W$5
	expxmx=exp(xbeta-mxbeta)

	// $X$5=EXP(U10-X2)
	expdxmdx=exp(dxbeta-mdxbeta)

	// building vector M - askla{i}==$W$5*(-LN(E2)) and askla{i}=F2*K2*$W$5*G2
	M1=expxmx*(-ln(E[1,1]))

	// F - acs**exp(x-mx)=POWER(E2,$W$5)
	F=E:^expxmx
	// K - ads**exp(dx-mdx)=POWER(J3,$X$5)
	K=J:^expdxmdx

	// M=askla{i}=F2*K2*$W$5*G2
	M2=expxmx*F[|.,.\rows(F)-1,.|]:*K[|.,.\rows(K)-1,.|]:*G
	M=(M1\M2)

	riskscore=colsum(M)
	riskvector[i]=riskscore

	//--------------------- Optimal risk score ---------------------

	// Q - make input vector into column vector
	Xopt=X[1]\X[2]\ln(110)\0\0\ln(22)\0

	// T26 - Xopt*beta
	xoptbeta=Xopt'*coeffx
	dxoptbeta=Xopt'*coeffdx

	// $W$26=EXP(T35-W2)
	expxoptmx=exp(xoptbeta-mxbeta)

	// $X$26=EXP(U35-X2)
	expdxoptmdx=exp(dxoptbeta-mdxbeta)

	// building vector AB - askla{1}_OPT=$W$26*(-LN(E2)) and askla{i}_OPT=Z2*AA2*$W$26*G2
	AB1=expxoptmx*(-ln(E[1,1]))
	// Z - acs**exp(x-mx) OPT=POWER(E2,$W$26)
	Z=E:^expxoptmx
	// AA - ads**exp(dx-mdx) OPT==POWER(J2,$X$26)
	AA=J:^expdxoptmdx
	// G - log{i-1}-log{i}, subset |top_left_row,top_left_column\bottom_right_row,bottom_right_column| - already made
	// AB=askla{i}_OPT=Z2*AA2*$W$26*G2
	AB2=expxoptmx*Z[|.,.\rows(Z)-1,.|]:*AA[|.,.\rows(AA)-1,.|]:*G
	AB=(AB1\AB2)
	riskscoreopt=colsum(AB)
	riskvectoropt[i]=riskscoreopt

	//--------------------- Normal risk score ---------------------

	// Q - make input vector into column vector
	Xnor=X[1]\X[2]\ln(125)\0\0\ln(22.5)\0

	// T26 - Xnor*beta
	xnorbeta=Xnor'*coeffx
	dxnorbeta=Xnor'*coeffdx

	// $W$39=EXP(T48-W2)
	expxnormx=exp(xnorbeta-mxbeta)

	// $X$39=EXP(U48-X2)
	expdxnormdx=exp(dxnorbeta-mdxbeta)

	// building vector AB - askla{1}_OPT=$W$26*(-LN(E2)) and askla{i}_OPT=Z2*AA2*$W$26*G2
	// building vector AF - askla{1}_NOR=$W$39*(-LN(E2)) and askla{i}_NOR=AD2*AE2*$W$39*G2
	AF1=expxnormx*(-ln(E[1,1]))
	// AD - acs**exp(x-mx) NOR=POWER(E2,$W$39)
	AD=E:^expxnormx
	// AE - ads**exp(dx-mdx) NOR=POWER(J2,$X$39)
	AE=J:^expdxnormdx
	// G - log{i-1}-log{i}, subset |top_left_row,top_left_column\bottom_right_row,bottom_right_column| - already made
	//askla{i}_NOR=AD2*AE2*$W$39*G2
	AF2=expxnormx*AD[|.,.\rows(AD)-1,.|]:*AE[|.,.\rows(AE)-1,.|]:*G
	AF=(AF1\AF2)
	riskscorenor=colsum(AF)
	riskvectornor[i]=riskscorenor

	}

// export risk score vectors into pre-created Stata variables

st_store(., "friskcvd30", touse, riskvector)
st_store(., "friskcvd30opt", touse, riskvectoropt)
st_store(., "friskcvd30nor", touse, riskvectornor)

}

end



//---------------------- Mata sub-program computing hard CVD risk ------------------------------

mata

real matrix RISKH (string scalar male, string scalar age, string scalar sbp, string scalar smoke, string scalar trtbp, string scalar bmi, string scalar diab, string scalar touse)
{
real scalar N, xbeta
real vector _sex, _logage, _logsbp, _smoke, _trtbp, _logbmi, _newdiab, _inputvars,  riskvector, riskvectoropt, riskvectornor

// create matrix which is a view on the current Stata dataset

st_view(inputvars=., ., (tokens(male), tokens(age), tokens(sbp), tokens(smoke), tokens(trtbp), tokens(bmi), tokens(diab)), touse)

// inputvars[1,.]

// Transform inputed values: sex, logage, logsbp, smoke, trtbp, logbmi, newdiab

_sex=inputvars[.,1]
_logage=ln(inputvars[.,2])
_logsbp=ln(inputvars[.,3])
_smoke=inputvars[.,4]
_trtbp=inputvars[.,5]
_logbmi=ln(inputvars[.,6])
_newdiab=inputvars[.,7]

_inputvars=_sex,_logage,_logsbp,_smoke,_trtbp,_logbmi,_newdiab

//_inputvars[1,.]

// Load initial data (it's loaded in pieces due to Stata memory limitation)

// Load data for hard CVD risk

#delimit;

E1=(0.999886630525445\0.999886630525445\0.999886630525445\0.999773181783610\0.999659707520039\0.999659707520039\0.999546126999023\0.999546126999023\0.999546126999023\0.999432392958266\0.999432392958266\0.999318617454709\0.999318617454709\0.999318617454709\0.999204686675588\0.999090701396040\0.998976654298056\0.998862521415169\0.998862521415169\0.998862521415169\0.998748246815426\0.998748246815426\0.998633851568511\0.998519450607380\0.998519450607380\0.998519450607380\0.998404804497814\0.998290128216189\0.998175403636370\0.998175403636370\
0.998060640541926\0.998060640541926\0.997945825344778\0.997945825344778\0.997945825344778\0.997945825344778\0.997830583829021\0.997715314550177\0.997600013880689\0.997600013880689\0.997484514067297\0.997368649957130\0.997252698726172\0.997136705915789\0.997020684655711\0.996904571879583\0.996788278084259\0.996671865766415\0.996671865766415\0.996671865766415\0.996555198095581\0.996438383096915\0.996321272457265\0.996321272457265\0.996321272457265\0.996203964097697\0.996086545698374\0.995969071148810\0.995851565094900\0.995734062149094\
0.995616497533746\0.995498816949754\0.995380943915757\0.995263001108616\0.995144774125882\0.995026507430955\0.994908251817238\0.994789805199607\0.994789805199607\0.994671351427472\0.994552877608217\0.994434352114371\0.994315831255141\0.994315831255141\0.994197209524056\0.994197209524056\0.994078259921501\0.994078259921501\0.994078259921501\0.994078259921501\0.993959039923569\0.993839802688272\0.993720550621687\0.993601280296658\0.993601280296658\0.993601280296658\0.993481808019070\0.993362157433778\0.993242397471890\0.993122593212210\
0.993002760978187\0.993002760978187\0.993002760978187\0.993002760978187\0.993002760978187\0.992882534522211\0.992762246020835\0.992641949555576\0.992521648119343\0.992401329111977\0.992281004676486\0.992281004676486\0.992160588972585\0.992160588972585\0.992160588972585\0.992039739862808\0.992039739862808\0.991918856784605\0.991797937369938\0.991676992638398\0.991676992638398\0.991555951025243\0.991434842543138\0.991313726577344\0.991192515315033\0.991071279246091\0.990950038622529\0.990828722661927\0.990707374996135\0.990585981616336\
0.990464373932848\0.990464373932848\0.990342715536429\0.990342715536429\0.990220902821429\0.990099047566470\0.989977186174996\0.989977186174996\0.989854921380006\0.989732647493911\0.989732647493911\0.989610110803402\0.989487472970173\0.989487472970173\0.989364778903279\0.989242079687464\0.989242079687464\0.989119143615312\0.988996138012341\0.988996138012341\0.988873102534730\0.988750048953147\0.988750048953147\0.988626870835478\0.988503637169224\0.988380391607727\0.988256950564465\0.988133448086608\0.988009943958225\0.988009943958225\
0.987885868746271\0.987761635589049\0.987761635589049\0.987637298893029\0.987512930365283\0.987512930365283\0.987388451082612\0.987263963884482\0.987139413257910\0.987139413257910\0.987014839595256\0.987014839595256\0.987014839595256\0.986890087466055\0.986765302916371\0.986640487947397\0.986640487947397\0.986640487947397\0.986515077488795\0.986389511195217\0.986389511195217\0.986263795209101\0.986263795209101\0.986137964727328\0.986012080240355\0.985886201323262\0.985760190409971\0.985760190409971\0.985634017315812\0.985634017315812\
0.985634017315812\0.985507597071419\0.985381063299762\0.985381063299762\0.985254334489131\0.985254334489131\0.985127480334990\0.985000514631362\0.984873438261428\0.984873438261428\0.984745895836353\0.984618336069284\0.984490686994074\0.984490686994074\0.984490686994074\0.984362302279141\0.984233898641617\0.984233898641617\0.984233898641617\0.984105093289864\0.984105093289864\0.984105093289864\0.983976102358288\0.983976102358288\0.983847023848751\0.983717927423257\0.983588774317126\0.983459603366769\0.983330384835341\0.983330384835341\
0.983200981065941\0.983200981065941\0.983200981065941\0.983071492015186\0.983071492015186\0.983071492015186\0.982941761459511\0.982941761459511\0.982811954173344\0.982682146487368\0.982552336511119\0.982422463206189\0.982292552045872\0.982162644533964\0.982162644533964\0.982032680187125\0.981902697900440\0.981772696977598\0.981772696977598\0.981642524915069\0.981642524915069\0.981512299921484\0.981381990712614\0.981381990712614\0.981381990712614\0.981251490016707\0.981251490016707\0.981120918888746\0.980990204070858\0.980859434884022);
E2=(0.980728595871431\0.980597722728521\0.980597722728521\0.980466721278427\0.980335693906936\0.980335693906936\0.980204585476906\0.980204585476906\0.980204585476906\0.980072912011481\0.980072912011481\0.979941063061451\0.979941063061451\0.979809124823225\0.979809124823225\0.979677038555460\0.979677038555460\0.979677038555460\0.979544517847708\0.979544517847708\0.979411754150162\0.979278909781664\0.979146068014201\0.979146068014201\0.979146068014201\0.979012936819995\0.979012936819995\0.978879628229522\0.978746282145256\0.978612884195268\
0.978612884195268\0.978479425629257\0.978345938089983\0.978345938089983\0.978345938089983\0.978345938089983\0.978212217853033\0.978212217853033\0.978212217853033\0.978212217853033\0.978078262013495\0.977944260955715\0.977810249428753\0.977676191474521\0.977541957529621\0.977541957529621\0.977407567322684\0.977407567322684\0.977273050441480\0.977273050441480\0.977138406868138\0.977003739323685\0.977003739323685\0.976868743344764\0.976868743344764\0.976868743344764\0.976733635471050\0.976733635471050\0.976733635471050\0.976598298286308\
0.976462869992546\0.976462869992546\0.976462869992546\0.976327021183305\0.976191176880041\0.976055221040328\0.975919226234551\0.975783081457982\0.975783081457982\0.975783081457982\0.975646848628817\0.975510467559389\0.975510467559389\0.975373686793582\0.975236869841534\0.975236869841534\0.975099948661104\0.975099948661104\0.975099948661104\0.975099948661103\0.974962505085787\0.974824556962316\0.974686407414055\0.974548240738906\0.974410042206478\0.974410042206478\0.974410042206478\0.974271709851502\0.974133299124362\0.973994778213898\
0.973994778213898\0.973856177637587\0.973856177637587\0.973717398024011\0.973717398024011\0.973578515901090\0.973439647667139\0.973439647667139\0.973439647667139\0.973300604710232\0.973300604710232\0.973161525227659\0.973022433614053\0.972883257331483\0.972883257331483\0.972744010802200\0.972604208661615\0.972604208661615\0.972464256190192\0.972464256190192\0.972324226265382\0.972184203074193\0.972044110580574\0.971903892610097\0.971763630419243\0.971623175913210\0.971623175913210\0.971482517680803\0.971341745014287\0.971200968149411\
0.971060168180826\0.971060168180825\0.970919097344585\0.970778007715525\0.970636925387611\0.970636925387611\0.970495752008433\0.970354575511489\0.970213344703287\0.970071850234189\0.969930191108740\0.969930191108740\0.969930191108740\0.969930191108740\0.969787849412495\0.969645484518641\0.969645484518641\0.969502565138095\0.969502565138095\0.969359523735936\0.969216469806275\0.969216469806275\0.969073342917074\0.969073342917074\0.968929949632928\0.968786516251905\0.968786516251905\0.968786516251905\0.968642927513320\0.968499275416237\
0.968499275416237\0.968355469579535\0.968355469579535\0.968211532154975\0.968067459931513\0.968067459931513\0.967923010835768\0.967778483961758\0.967633837743092\0.967489096119651\0.967344294985510\0.967344294985510\0.967199305025095\0.967054247380810\0.967054247380810\0.967054247380810\0.967054247380810\0.966909023336799\0.966763748200006\0.966763748200006\0.966618372034198\0.966472929517210\0.966327498577099\0.966181852619510\0.966181852619510\0.966181852619510\0.966181852619510\0.966181852619510\0.966035945887102\0.966035945887102\
0.966035945887102\0.966035945887102\0.966035945887102\0.965889721494597\0.965889721494597\0.965889721494597\0.965743192637638\0.965743192637638\0.965596537819963\0.965596537819963\0.965596537819963\0.965449652473018\0.965302661265722\0.965155672902541\0.965155672902541\0.965155672902541\0.965008347769582\0.965008347769582\0.964860987233673\0.964713570142889\0.964566104598991\0.964418254800283\0.964418254800283\0.964270229317861\0.964270229317861\0.964121776103402\0.963973267240205\0.963973267240205\0.963824707681306\0.963676126151817\
0.963527337676703\0.963527337676703\0.963527337676703\0.963527337676703\0.963378317421182\0.963378317421182\0.963378317421182\0.963378317421182\0.963378317421182\0.963228013037217\0.963228013037217\0.963228013037217\0.963228013037216\0.963228013037216\0.963077194214014\0.962926222269167\0.962775245066506\0.962775245066506\0.962624134289052\0.962472921585934\0.962472921585934\0.962321576095935\0.962170232232869\0.962018904842703\0.962018904842703\0.961867483171417\0.961715912495578\0.961564113283023\0.961412277115300\0.961260102837572);
E3=(0.961260102837572\0.961107718616193\0.960955306355814\0.960802878260876\0.960650323571821\0.960650323571821\0.960497272253471\0.960344200864374\0.960344200864374\0.960190819903857\0.960190819903856\0.960037356561347\0.959883868438120\0.959730375829100\0.959576293313762\0.959422204634417\0.959422204634417\0.959422204634417\0.959268012713151\0.959268012713151\0.959113818447173\0.958959574767193\0.958805307286350\0.958805307286350\0.958805307286350\0.958805307286350\0.958805307286350\0.958805307286350\0.958805307286350\0.958649935771748\
0.958494530895195\0.958338879257491\0.958338879257491\0.958338879257491\0.958338879257491\0.958338879257491\0.958182476090919\0.958026046036376\0.958026046036376\0.958026046036376\0.958026046036376\0.958026046036376\0.958026046036376\0.958026046036376\0.958026046036376\0.957868817695577\0.957868817695577\0.957710726672411\0.957710726672411\0.957552480931078\0.957394231680967\0.957235892682486\0.957077569025085\0.957077569025085\0.957077569025085\0.957077569025085\0.956919028230479\0.956760398598974\0.956760398598974\0.956601719154726\
0.956601719154726\0.956601719154726\0.956601719154726\0.956601719154726\0.956442722255304\0.956442722255304\0.956442722255304\0.956283449511584\0.956123916528599\0.956123916528599\0.956123916528599\0.956123916528599\0.956123916528599\0.955964064535591\0.955964064535591\0.955964064535591\0.955804076879063\0.955804076879063\0.955804076879063\0.955643679748075\0.955643679748075\0.955483107752460\0.955483107752460\0.955483107752460\0.955322359161902\0.955322359161902\0.955161581870172\0.955000797596301\0.955000797596301\0.955000797596301\
0.954839879904012\0.954678743840369\0.954678743840369\0.954517576557785\0.954356364590645\0.954195104288052\0.954033679441184\0.953872234892689\0.953710721452756\0.953548984184971\0.953548984184971\0.953387126036164\0.953387126036164\0.953387126036164\0.953387126036164\0.953387126036164\0.953387126036164\0.953223895890583\0.953060098213690\0.952896185398308\0.952896185398308\0.952896185398308\0.952896185398308\0.952732196682106\0.952568124387221\0.952568124387221\0.952568124387221\0.952568124387221\0.952403966515807\0.952403966515807\
0.952239790267357\0.952239790267357\0.952075392737739\0.952075392737739\0.951910724271994\0.951746039638923\0.951581365774534\0.951581365774534\0.951416534251524\0.951416534251524\0.951416534251524\0.951251445776165\0.951086341203278\0.950921236267035\0.950921236267035\0.950755882321241\0.950590372862574\0.950590372862574\0.950590372862574\0.950590372862574\0.950424544299311\0.950258700234950\0.950258700234950\0.950092453744028\0.950092453744028\0.949926185416667\0.949926185416667\0.949759844937253\0.949593459047046\0.949427028642141\
0.949260571338277\0.949260571338277\0.949094060270738\0.948927457078052\0.948760844216889\0.948760844216889\0.948760844216889\0.948593699546575\0.948593699546575\0.948593699546575\0.948426195233889\0.948258686945681\0.948090147671575\0.947921452166076\0.947752748069107\0.947752748069107\0.947752748069107\0.947752748069107\0.947583265817199\0.947583265817199\0.947583265817199\0.947413642314263\0.947244000360976\0.947074234153705\0.947074234153705\0.946904382888881\0.946904382888881\0.946904382888881\0.946734379286866\0.946734379286866\
0.946734379286866\0.946734379286866\0.946564223613675\0.946564223613675\0.946393328124295\0.946222339686086\0.946222339686086\0.946222339686086\0.946051103971502\0.946051103971502\0.946051103971502\0.945879648126614\0.945708191632689\0.945708191632689\0.945536623156267\0.945365058152364\0.945365058152364\0.945365058152364\0.945193265225276\0.945021368120908\0.944849407674976\0.944849407674976\0.944849407674976\0.944677295751086\0.944505150762416\0.944332944595276\0.944160749127808\0.944160749127808\0.944160749127808\0.943988321857951\
0.943815787649749\0.943815787649749\0.943643103134062\0.943470373716900\0.943297567174678\0.943124697961965\0.943124697961965\0.943124697961965\0.942951442107339\0.942951442107339\0.942777851061869\0.942777851061869\0.942604052016978\0.942604052016978\0.942430067047293\0.942256058092436\0.942256058092436\0.942256058092436\0.942256058092436\0.942256058092436\0.942256058092436\0.942081356199637\0.941906395712068\0.941906395712068\0.941731246902090\0.941731246902089\0.941556068099271\0.941380854380402\0.941380854380402\0.941205472352831);
E4=(0.941205472352831\0.941029943818062\0.940854440903021\0.940854440903021\0.940678700139399\0.940678700139399\0.940678700139399\0.940678700139399\0.940502645047277\0.940502645047277\0.940502645047277\0.940326317095853\0.940149936989458\0.939973453516870\0.939796956041275\0.939796956041275\0.939620305100016\0.939620305100016\0.939443420580542\0.939266539291356\0.939089644600921\0.938912766759595\0.938912766759595\0.938735463633331\0.938735463633331\0.938557783358138\0.938379854964084\0.938379854964084\0.938379854964084\0.938379854964084\
0.938379854964084\0.938379854964084\0.938201374720575\0.938022885626913\0.937844403820960\0.937844403820960\0.937665694787405\0.937665694787405\0.937665694787405\0.937486619757432\0.937307441432943\0.937307441432943\0.937307441432943\0.937307441432943\0.937307441432943\0.937307441432943\0.937307441432943\0.937307441432943\0.937127527351502\0.937127527351502\0.937127527351502\0.936947454078083\0.936947454078083\0.936767350046493\0.936767350046493\0.936587001727456\0.936587001727456\0.936587001727456\0.936587001727456\0.936406548298852\
0.936406548298852\0.936406548298852\0.936225759990117\0.936045000365154\0.936045000365154\0.935864151110763\0.935683316928326\0.935683316928326\0.935683316928326\0.935502062987408\0.935320700583546\0.935320700583546\0.935139160740455\0.934957288992079\0.934775350749166\0.934775350749166\0.934775350749166\0.934593269359590\0.934593269359590\0.934411104701699\0.934228973443294\0.934228973443294\0.934228973443294\0.934228973443294\0.934228973443294\0.934045692327439\0.934045692327439\0.934045692327439\0.933862298847492\0.933862298847492\
0.933678218412528\0.933678218412528\0.933678218412528\0.933678218412528\0.933678218412528\0.933493512835703\0.933493512835703\0.933493512835703\0.933308304849000\0.933123068785316\0.932937764274492\0.932937764274492\0.932937764274492\0.932937764274492\0.932937764274492\0.932752079533806\0.932752079533806\0.932752079533806\0.932752079533806\0.932565600577627\0.932565600577627\0.932379032441629\0.932379032441629\0.932192415546694\0.932192415546694\0.932005714215769\0.932005714215770\0.931818998032429\0.931818998032429\0.931818998032429\
0.931818998032429\0.931631429936796\0.931631429936796\0.931443628389004\0.931255750735270\0.931255750735270\0.931067543664211\0.931067543664211\0.930879319117553\0.930690789718208\0.930690789718209\0.930501609098186\0.930312406273162\0.930123171250006\0.929933857739539\0.929933857739539\0.929933857739539\0.929933857739539\0.929743733344613\0.929553097852639\0.929553097852639\0.929362407163661\0.929171650648021\0.928980882097039\0.928789999147378\0.928789999147378\0.928789999147378\0.928789999147378\0.928789999147378\0.928598604231665\
0.928407193773450\0.928407193773450\0.928407193773450\0.928215033201763\0.928215033201763\0.928215033201763\0.928022401809240\0.927829578617774\0.927829578617774\0.927829578617774\0.927829578617774\0.927636213601494\0.927442758397469\0.927442758397469\0.927442758397469\0.927442758397469\0.927248290356659\0.927053731725893\0.927053731725893\0.927053731725893\0.926858747795846\0.926858747795846\0.926858747795846\0.926858747795846\0.926858747795846\0.926858747795846\0.926661889178441\0.926465051596524\0.926465051596524\0.926465051596524\
0.926268058215580\0.926070642574867\0.926070642574868\0.926070642574867\0.926070642574868\0.925872581408546\0.925674548125803\0.925674548125803\0.925674548125803\0.925674548125803\0.925674548125803\0.925674548125803\0.925674548125803\0.925674548125803\0.925674548125803\0.925674548125803\0.925674548125803\0.925674548125803\0.925474528262458\0.925474528262458\0.925474528262458\0.925474528262458\0.925274197036722\0.925274197036722\0.925073608084967\0.924872954147309\0.924872954147309\0.924872954147309\0.924872954147309\0.924671334079094\
0.924671334079094\0.924671334079094\0.924671334079094\0.924469210951846\0.924265952899392\0.924062577731620\0.924062577731620\0.923859098893356\0.923655555774034\0.923655555774034\0.923451994931359\0.923451994931359\0.923247945406606\0.923043841389494\0.922839309284649\0.922634733517995\0.922430132698179\0.922430132698179\0.922225300255632\0.922020200729106\0.921815078107734\0.921609733327381\0.921609733327381\0.921609733327381\0.921404111935558\0.921404111935558\0.921404111935558\0.921404111935558\0.921404111935558\0.921198024397284);
E5=(0.920991786843172\0.920785378149555\0.920785378149555\0.920785378149555\0.920785378149555\0.920578149688686\0.920578149688686\0.920370451354702\0.920370451354703\0.920370451354703\0.920161653056505\0.919952816845927\0.919952816845927\0.919743759718265\0.919743759718265\0.919743759718265\0.919533671140828\0.919323591278294\0.919323591278294\0.919323591278294\0.919112785910033\0.919112785910033\0.919112785910033\0.918901744960664\0.918690677754476\0.918479598557798\0.918479598557798\0.918479598557798\0.918479598557798\0.918268031311692\
0.918268031311692\0.918268031311692\0.918056256118922\0.917844424717364\0.917632194111213\0.917419977652232\0.917419977652232\0.917207312486977\0.916994319103474\0.916781278484560\0.916781278484560\0.916568163958661\0.916355076631379\0.916355076631379\0.916141835003065\0.916141835003065\0.915928368754571\0.915714779538503\0.915500794700247\0.915286373338954\0.915286373338954\0.915286373338954\0.915071573325881\0.915071573325881\0.914856511565752\0.914856511565752\0.914856511565752\0.914856511565752\0.914640701200380\0.914640701200380\
0.914640701200380\0.914640701200379\0.914424260074668\0.914424260074668\0.914424260074668\0.914206902756272\0.913988811150943\0.913769894159457\0.913769894159457\0.913769894159457\0.913769894159457\0.913549463626504\0.913549463626504\0.913549463626504\0.913328394917753\0.913107356063671\0.913107356063671\0.913107356063671\0.913107356063671\0.912884333974300\0.912884333974300\0.912884333974300\0.912884333974300\0.912660192079330\0.912660192079330\0.912435468279400\0.912435468279400\0.912209995209672\0.912209995209673\0.911984033420350\
0.911984033420350\0.911984033420350\0.911984033420350);

J1=(1.000000000000000\0.999877382315385\0.999754739237093\0.999754739237093\0.999754739237093\0.999631956390161\0.999631956390161\0.999509066575281\0.999386079705851\0.999386079705851\0.999263006872402\0.999263006872402\0.999139888977779\0.999016678346870\0.999016678346870\0.999016678346870\0.999016678346870\0.999016678346870\0.998892976437443\0.998769283598168\0.998769283598168\0.998645434338833\0.998645434338833\0.998645434338833\0.998521445699783\0.998397292063546\0.998397292063546\0.998397292063546\0.998397292063546\0.998272944972888\
0.998272944972888\0.998148505546113\0.998148505546113\0.998023964637041\0.997899177258379\0.997774397828149\0.997774397828149\0.997774397828149\0.997774397828149\0.997649329665411\0.997649329665411\0.997649329665411\0.997649329665411\0.997649329665411\0.997649329665411\0.997649329665411\0.997649329665411\0.997649329665411\0.997523172716857\0.997397001753030\0.997397001753030\0.997397001753030\0.997397001753030\0.997270314447709\0.997143622212036\0.997143622212036\0.997143622212036\0.997143622212036\0.997143622212036\0.997143622212036\
0.997143622212036\0.997143622212036\0.997143622212036\0.997143622212036\0.997143622212036\0.997143622212036\0.997143622212036\0.997143622212036\0.997015747936520\0.997015747936520\0.997015747936520\0.997015747936520\0.997015747936520\0.996887677957909\0.996887677957909\0.996759540986866\0.996759540986866\0.996631054277965\0.996502453166396\0.996373863924988\0.996373863924988\0.996373863924988\0.996373863924988\0.996373863924988\0.996245105788307\0.996116315814913\0.996116315814913\0.996116315814913\0.996116315814913\0.996116315814913\
0.996116315814913\0.995986839734619\0.995857361145062\0.995727823821760\0.995598175589944\0.995598175589944\0.995598175589944\0.995598175589944\0.995598175589944\0.995598175589944\0.995598175589944\0.995468184047067\0.995468184047067\0.995337899800486\0.995207582891355\0.995207582891355\0.995077140357245\0.995077140357245\0.995077140357245\0.995077140357245\0.994946542640124\0.994946542640124\0.994946542640124\0.994946542640124\0.994946542640124\0.994946542640124\0.994946542640124\0.994946542640124\0.994946542640124\0.994946542640124\
0.994946542640124\0.994815157768099\0.994815157768098\0.994683627527736\0.994683627527736\0.994683627527736\0.994683627527736\0.994551531000362\0.994551531000362\0.994551531000362\0.994419352352335\0.994419352352335\0.994419352352335\0.994286788753509\0.994286788753509\0.994286788753509\0.994153920135270\0.994153920135270\0.994153920135270\0.994020860490313\0.994020860490313\0.994020860490313\0.993887622820355\0.993887622820355\0.993887622820355\0.993887622820355\0.993887622820355\0.993887622820355\0.993887622820355\0.993753885903036\
0.993753885903036\0.993753885903036\0.993619608697495\0.993619608697495\0.993619608697495\0.993485251583580\0.993485251583580\0.993485251583580\0.993485251583580\0.993350708253490\0.993350708253490\0.993216101218718\0.993081328382637\0.993081328382637\0.993081328382637\0.993081328382637\0.992946352547719\0.992810951301085\0.992810951301085\0.992810951301085\0.992675218495593\0.992675218495593\0.992539220327278\0.992539220327278\0.992539220327278\0.992539220327278\0.992539220327278\0.992402997645354\0.992402997645354\0.992266652588372\
0.992130258730334\0.992130258730334\0.992130258730334\0.991993578257334\0.991993578257334\0.991856716918194\0.991856716918194\0.991856716918194\0.991856716918194\0.991719468854811\0.991719468854811\0.991719468854811\0.991719468854811\0.991581685315082\0.991443628156897\0.991443628156897\0.991443628156897\0.991305379654443\0.991167072429714\0.991167072429714\0.991028617071443\0.990890175441504\0.990890175441504\0.990751558100540\0.990751558100540\0.990751558100540\0.990751558100540\0.990751558100540\0.990751558100540\0.990612634352295\
0.990612634352295\0.990473621305981\0.990334593938515\0.990334593938515\0.990195460016739\0.990056334787741\0.990056334787741\0.989917027423623\0.989917027423623\0.989917027423623\0.989917027423623\0.989917027423623\0.989917027423623\0.989917027423623\0.989777425967349\0.989777425967349\0.989777425967349\0.989777425967349\0.989637645335595\0.989637645335595\0.989497767031296\0.989497767031296\0.989497767031296\0.989357669551902\0.989217477891748\0.989217477891748\0.989077219114354\0.989077219114354\0.989077219114354\0.989077219114354);
J2=(0.989077219114354\0.989077219114354\0.988936512821478\0.988936512821478\0.988936512821478\0.988795593078810\0.988795593078810\0.988654524331113\0.988513437257583\0.988513437257583\0.988372072871561\0.988372072871561\0.988230595438959\0.988230595438959\0.988088949012826\0.988088949012826\0.987946787984170\0.987804572124083\0.987804572124083\0.987662108496535\0.987662108496535\0.987662108496535\0.987662108496535\0.987519331076970\0.987376548328310\0.987376548328310\0.987233569574694\0.987233569574694\0.987233569574694\0.987233569574694\
0.987090264656929\0.987090264656929\0.987090264656929\0.986946732860676\0.986803123182412\0.986659485688719\0.986659485688719\0.986515671636055\0.986371859100067\0.986227962156667\0.986227962156667\0.986227962156667\0.986227962156667\0.986227962156667\0.986227962156667\0.986083684883257\0.986083684883257\0.985939256131526\0.985939256131526\0.985794626894039\0.985794626894039\0.985794626894039\0.985649735456292\0.985649735456292\0.985504739124788\0.985359713958848\0.985359713958848\0.985214547114021\0.985069356411770\0.985069356411770\
0.985069356411770\0.984923573431527\0.984777568349374\0.984777568349374\0.984777568349374\0.984777568349374\0.984777568349374\0.984777568349374\0.984631210546111\0.984484815787104\0.984484815787104\0.984484815787104\0.984338113813830\0.984338113813830\0.984338113813830\0.984191175992438\0.984191175992438\0.984044011984509\0.983896794135830\0.983749534097872\0.983749534097872\0.983749534097872\0.983749534097872\0.983749534097872\0.983749534097872\0.983601555651326\0.983453568321210\0.983453568321210\0.983453568321210\0.983453568321210\
0.983305154273215\0.983305154273214\0.983156628490682\0.983156628490682\0.983007969216944\0.983007969216944\0.983007969216944\0.982859200025887\0.982710424503402\0.982710424503402\0.982561566754414\0.982561566754414\0.982561566754414\0.982561566754414\0.982412557657166\0.982412557657166\0.982412557657166\0.982263267586500\0.982263267586500\0.982113763014034\0.982113763014034\0.982113763014034\0.982113763014034\0.982113763014034\0.982113763014034\0.982113763014034\0.981963533645999\0.981963533645999\0.981963533645999\0.981963533645999\
0.981963533645999\0.981812906045177\0.981812906045177\0.981812906045177\0.981812906045177\0.981661945147075\0.981661945147075\0.981661945147075\0.981661945147075\0.981661945147075\0.981661945147075\0.981510210581296\0.981358377980760\0.981206512932880\0.981206512932880\0.981206512932880\0.981054374634427\0.981054374634427\0.980901912123267\0.980901912123267\0.980901912123267\0.980749267026135\0.980749267026135\0.980596438004534\0.980596438004534\0.980596438004534\0.980443421687700\0.980290408028649\0.980290408028649\0.980290408028649\
0.980137157129925\0.980137157129925\0.979983794024871\0.979983794024871\0.979983794024871\0.979830162904288\0.979830162904288\0.979830162904288\0.979830162904288\0.979830162904288\0.979830162904288\0.979675974030892\0.979675974030892\0.979675974030892\0.979521382098135\0.979366801061457\0.979212224565338\0.979212224565338\0.979212224565338\0.979057490865897\0.979057490865897\0.979057490865897\0.979057490865897\0.979057490865897\0.978902556870616\0.978747593362878\0.978592558989321\0.978437478917565\0.978437478917565\0.978282209022862\
0.978126898404028\0.977971592225357\0.977816211215451\0.977816211215451\0.977660465461769\0.977504634881541\0.977504634881542\0.977348616363237\0.977348616363237\0.977192543126441\0.977036362414105\0.977036362414105\0.977036362414105\0.977036362414105\0.976880014467439\0.976723666965032\0.976723666965032\0.976567074907804\0.976567074907804\0.976567074907804\0.976567074907804\0.976567074907804\0.976410028823628\0.976410028823628\0.976252669978633\0.976252669978633\0.976252669978633\0.976095105712386\0.976095105712386\0.976095105712386\
0.976095105712386\0.975937294153744\0.975779279095949\0.975621258061589\0.975621258061589\0.975462583381072\0.975303688695002\0.975144647178305\0.974985340121587\0.974985340121587\0.974825835177759\0.974666275961057\0.974506488659083\0.974346697605226\0.974346697605226\0.974346697605226\0.974346697605226\0.974186334789237\0.974186334789237\0.974186334789237\0.974025828179885\0.974025828179884\0.974025828179885\0.974025828179885\0.973865136889873\0.973865136889873\0.973865136889873\0.973865136889873\0.973865136889873\0.973865136889873);
J3=(0.973703708178830\0.973703708178830\0.973703708178830\0.973703708178830\0.973703708178830\0.973541461350011\0.973541461350011\0.973541461350011\0.973378703717403\0.973378703717403\0.973215826673056\0.973215826673056\0.973215826673056\0.973215826673056\0.973215826673056\0.973215826673056\0.973052417443192\0.972888992210845\0.972888992210845\0.972725547073634\0.972725547073634\0.972725547073634\0.972725547073634\0.972561904326310\0.972398143893568\0.972234165755981\0.972070187198272\0.971906060540717\0.971741794319356\0.971741794319356\
0.971741794319356\0.971741794319356\0.971576917939075\0.971411807392104\0.971246446288561\0.971081014742375\0.971081014742375\0.971081014742375\0.970915377665463\0.970749650861005\0.970583889109737\0.970418117398808\0.970252289428832\0.970086210928904\0.969920145771451\0.969920145771451\0.969753574894551\0.969753574894551\0.969586401950765\0.969586401950765\0.969586401950765\0.969586401950765\0.969586401950765\0.969418925165711\0.969251466328004\0.969083980285280\0.969083980285280\0.969083980285280\0.968916317991869\0.968916317991869\
0.968748598121546\0.968580878549095\0.968413050577942\0.968245214941023\0.968245214941023\0.968077236623766\0.967909131802710\0.967909131802710\0.967909131802710\0.967740501251592\0.967571770061472\0.967403013274084\0.967234239958707\0.967234239958707\0.967065308864984\0.966896385303187\0.966896385303187\0.966727349294869\0.966558317663835\0.966558317663835\0.966389028769427\0.966389028769427\0.966219573318300\0.966049998856747\0.966049998856747\0.965880379509070\0.965880379509070\0.965880379509070\0.965710661169920\0.965540946035333\
0.965540946035333\0.965540946035333\0.965370754299050\0.965370754299050\0.965370754299050\0.965370754299050\0.965370754299050\0.965370754299050\0.965370754299050\0.965370754299050\0.965199787527636\0.965199787527636\0.965028704601583\0.964857287746553\0.964685837323690\0.964514161608962\0.964342506956207\0.964342506956207\0.964342506956207\0.964342506956207\0.964169531896089\0.963996543970675\0.963823579828715\0.963823579828715\0.963823579828715\0.963650468841604\0.963477361875273\0.963304235336023\0.963304235336023\0.963131096465925\
0.963131096465925\0.962957837976841\0.962957837976841\0.962784322031305\0.962784322031305\0.962784322031305\0.962784322031305\0.962610473729852\0.962610473729852\0.962436279832276\0.962262094951067\0.962262094951068\0.962262094951067\0.962262094951068\0.962087751914124\0.962087751914124\0.962087751914124\0.961912720971110\0.961737613619076\0.961562525267156\0.961562525267156\0.961562525267156\0.961387132381365\0.961387132381365\0.961211588395338\0.961211588395338\0.961035985767719\0.961035985767719\0.961035985767719\0.961035985767719\
0.961035985767719\0.960860047666178\0.960860047666178\0.960860047666178\0.960860047666178\0.960683955435401\0.960507608916188\0.960507608916188\0.960331133695281\0.960154596520646\0.960154596520646\0.960154596520646\0.960154596520646\0.960154596520646\0.960154596520646\0.959976938616596\0.959798958319262\0.959620979825081\0.959620979825081\0.959442986673994\0.959265012212499\0.959265012212499\0.959265012212499\0.959265012212499\0.959086630116221\0.959086630116221\0.958908204050802\0.958729787947037\0.958729787947037\0.958551220522066\
0.958372633670360\0.958194073970496\0.958194073970496\0.958014995463407\0.958014995463407\0.958014995463407\0.957835556134349\0.957656061385763\0.957656061385763\0.957476320718517\0.957296560655553\0.957296560655553\0.957296560655553\0.957116684245007\0.957116684245007\0.957116684245007\0.956936670698443\0.956756454300902\0.956756454300902\0.956756454300902\0.956756454300902\0.956575932521137\0.956395352215026\0.956395352215026\0.956395352215026\0.956395352215026\0.956395352215026\0.956214377575624\0.956033371500509\0.956033371500509\
0.956033371500508\0.955852143473998\0.955852143473998\0.955852143473998\0.955852143473998\0.955852143473998\0.955670520547735\0.955488839517210\0.955488839517210\0.955306878498552\0.955306878498552\0.955124623986018\0.955124623986018\0.954942189295016\0.954942189295016\0.954942189295016\0.954759441980690\0.954576497709827\0.954393375214054\0.954210240519152\0.954027083471315\0.954027083471315\0.954027083471315\0.953843391732980\0.953843391732980\0.953659610250826\0.953659610250826\0.953659610250826\0.953475582340610\0.953475582340610);
J4=(0.953291402198746\0.953291402198746\0.953291402198746\0.953107107574818\0.953107107574818\0.952922626607897\0.952738089941614\0.952553256685820\0.952553256685820\0.952368289731231\0.952183132787505\0.952183132787504\0.952183132787505\0.952183132787505\0.952183132787505\0.951997370749014\0.951997370749014\0.951811479942043\0.951811479942043\0.951811479942043\0.951811479942043\0.951811479942043\0.951625254953933\0.951625254953933\0.951438478941965\0.951438478941965\0.951438478941965\0.951251384482227\0.951064293574907\0.950877058935636\
0.950689809172153\0.950502512789198\0.950502512789198\0.950502512789198\0.950502512789198\0.950315017385709\0.950315017385709\0.950127319653130\0.949939632197418\0.949939632197418\0.949939632197418\0.949751294248673\0.949562853462727\0.949374442435989\0.949185944168219\0.948997391353092\0.948808754341884\0.948620120305231\0.948620120305231\0.948431418628924\0.948242743469254\0.948242743469254\0.948053920065868\0.948053920065868\0.947864986593081\0.947864986593081\0.947675912801324\0.947486792241033\0.947297682563024\0.947297682563024\
0.947108550210175\0.946919203749106\0.946919203749106\0.946919203749106\0.946729640620813\0.946729640620813\0.946729640620813\0.946539953361653\0.946350168150479\0.946350168150479\0.946350168150479\0.946160021135842\0.946160021135842\0.946160021135842\0.946160021135842\0.945969464518774\0.945778758630404\0.945778758630404\0.945587938598075\0.945587938598075\0.945587938598075\0.945396871338807\0.945205731700962\0.945014496721482\0.944823153856457\0.944823153856457\0.944631459367357\0.944439697765413\0.944439697765413\0.944247843832257\
0.944247843832257\0.944055480986861\0.943862998703275\0.943670474049101\0.943477803525064\0.943477803525064\0.943284979127113\0.943091860998315\0.943091860998315\0.943091860998315\0.943091860998315\0.942898216830883\0.942704576652458\0.942510846400247\0.942317026739292\0.942317026739292\0.942123018616505\0.941928618052765\0.941734110935843\0.941734110935843\0.941539388783547\0.941539388783547\0.941344551176867\0.941344551176867\0.941149720553436\0.941149720553436\0.940954792353922\0.940954792353922\0.940759558616909\0.940564088000740\
0.940368394491591\0.940368394491591\0.940172666191057\0.940172666191057\0.940172666191057\0.939976366421197\0.939976366421197\0.939779890502935\0.939779890502935\0.939779890502935\0.939582687262359\0.939582687262359\0.939582687262359\0.939582687262359\0.939582687262359\0.939385093741012\0.939187240068946\0.938989177973288\0.938989177973289\0.938989177973288\0.938790267130564\0.938790267130564\0.938790267130564\0.938790267130564\0.938790267130564\0.938591055199027\0.938391847456297\0.938192380832533\0.937992886522510\0.937992886522510\
0.937992886522510\0.937792670517370\0.937592468995621\0.937592468995621\0.937391871509223\0.937190866176701\0.937190866176701\0.937190866176701\0.936989394402064\0.936787951253132\0.936586519000650\0.936586519000650\0.936586519000650\0.936384794239185\0.936183085980087\0.935981283067797\0.935981283067797\0.935981283067797\0.935778796899950\0.935576324814780\0.935576324814780\0.935373377095168\0.935170257448275\0.934966627711866\0.934763005998565\0.934559288808697\0.934559288808697\0.934559288808697\0.934355410666419\0.934151499901847\
0.934151499901847\0.934151499901847\0.933946876095152\0.933742151483729\0.933537131567790\0.933537131567790\0.933537131567790\0.933332004998998\0.933126885763364\0.932921660509063\0.932716458047412\0.932510876824801\0.932305259105415\0.932099345362601\0.931893348639551\0.931687380177541\0.931481432737922\0.931275275311487\0.931275275311487\0.931069075051751\0.930862820495296\0.930656475105424\0.930656475105424\0.930449998501725\0.930449998501725\0.930449998501725\0.930243311552208\0.930036488810083\0.929829453792992\0.929829453792992\
0.929622208346582\0.929414946722466\0.929207530326662\0.929207530326662\0.929207530326662\0.929207530326662\0.928999257499280\0.928999257499280\0.928999257499280\0.928790842176191\0.928790842176191\0.928582022945592\0.928582022945592\0.928582022945592\0.928582022945592\0.928582022945592\0.928582022945592\0.928372046129379\0.928372046129378\0.928372046129379\0.928372046129379\0.928372046129379\0.928161369312829\0.927950597331561\0.927950597331561\0.927739748321159\0.927528893624456\0.927317861074674\0.927106618756137\0.927106618756137);
J5=(0.927106618756137\0.927106618756137\0.926894579770270\0.926682280263304\0.926469947589908\0.926469947589908\0.926257284155597\0.926257284155597\0.926044010559335\0.925830260837497\0.925830260837497\0.925830260837497\0.925616177495272\0.925616177495272\0.925401463901501\0.925186617471379\0.925186617471379\0.925186617471379\0.924971129836086\0.924755557874795\0.924755557874796\0.924539610977740\0.924323638051896\0.924323638051896\0.924323638051896\0.924323638051896\0.924107423516060\0.923891219640908\0.923674943850368\0.923674943850368\
0.923458500718814\0.923241926837872\0.923241926837872\0.923241926837872\0.923241926837872\0.923241926837872\0.923024827952549\0.923024827952549\0.923024827952549\0.923024827952549\0.922807066861285\0.922807066861285\0.922807066861285\0.922589139571138\0.922589139571138\0.922370893821665\0.922370893821665\0.922370893821665\0.922370893821665\0.922370893821665\0.922151616938659\0.921932103128272\0.921932103128272\0.921712340510784\0.921712340510784\0.921492116270046\0.921271712987450\0.921051158798440\0.921051158798440\0.920830407364890\
0.920609705732859\0.920388690016989\0.920388690016989\0.920167212847615\0.919945043479497\0.919945043479497\0.919945043479497\0.919945043479497\0.919720658085710\0.919495991727049\0.919271143622711\0.919271143622711\0.919045503254822\0.918819741716551\0.918819741716551\0.918819741716551\0.918593650774855\0.918367206846034\0.918140177205785\0.918140177205785\0.917912453859927\0.917683706270257\0.917454902131317\0.917454902131317\0.917225642354725\0.917225642354725\0.916995764858996\0.916995764858996\0.916765558184566\0.916765558184566\
0.916534760403720\0.916303773010106\0.916072611845584);

#delimit cr

// assemble pieces of initial data into single vectors

E=E1\E2\E3\E4\E5
J=J1\J2\J3\J4\J5

// Create matrix G (used below) - log{i-1}-log{i}, subset |top_left_row,top_left_column\bottom_right_row,bottom_right_column|
E1=E[|1,.\rows(E)-1,.|]
E2=E[|2,.\rows(E),.|]
G=ln(E1)-ln(E2)

// input initial values of parameters

coeffx=(0.73413\3.07187\1.79939\0.79437\0.3926\1.1166\1.03634)
coeffdx=(0.47275\3.55871\1.5204\0.96516\0.11331\-0.27572\0.45868)
mxbeta=23.9352581083579
mdxbeta=19.7973162240596

// Define new vector containing risk scores
riskvector=J(rows(_inputvars),1,.)
riskvectoropt=J(rows(_inputvars),1,.)
riskvectornor=J(rows(_inputvars),1,.)

// loop over observations and compute risk score for each one

for (i=1; i<=rows(_inputvars); i++) {

	// select observation i of input data
	Xrow=_inputvars[| i,1\i,.|]

	//--------------------- Risk score ---------------------

	// Q - make input vector into column vector
	X=Xrow'

	// R - X*beta
	xbeta=X'*coeffx
	dxbeta=X'*coeffdx

	// $W$5
	expxmx=exp(xbeta-mxbeta)

	// $X$5=EXP(U10-X2)
	expdxmdx=exp(dxbeta-mdxbeta)

	// building vector M - askla{i}==$W$5*(-LN(E2)) and askla{i}=F2*K2*$W$5*G2
	M1=expxmx*(-ln(E[1,1]))

	// F - acs**exp(x-mx)=POWER(E2,$W$5)
	F=E:^expxmx
	// K - ads**exp(dx-mdx)=POWER(J3,$X$5)
	K=J:^expdxmdx

	// M=askla{i}=F2*K2*$W$5*G2
	M2=expxmx*F[|.,.\rows(F)-1,.|]:*K[|.,.\rows(K)-1,.|]:*G
	M=(M1\M2)
	riskscore=colsum(M)
	riskvector[i]=riskscore

	//--------------------- Optimal risk score ---------------------

	// Q - make input vector into column vector
	Xopt=X[1]\X[2]\ln(110)\0\0\ln(22)\0

	// T26 - Xopt*beta
	xoptbeta=Xopt'*coeffx
	dxoptbeta=Xopt'*coeffdx

	// $W$26=EXP(T35-W2)
	expxoptmx=exp(xoptbeta-mxbeta)

	// $X$26=EXP(U35-X2)
	expdxoptmdx=exp(dxoptbeta-mdxbeta)

	// building vector AB - askla{1}_OPT=$W$26*(-LN(E2)) and askla{i}_OPT=Z2*AA2*$W$26*G2
	AB1=expxoptmx*(-ln(E[1,1]))
	// Z - acs**exp(x-mx) OPT=POWER(E2,$W$26)
	Z=E:^expxoptmx
	// AA - ads**exp(dx-mdx) OPT==POWER(J2,$X$26)
	AA=J:^expdxoptmdx
	// G - log{i-1}-log{i}, subset |top_left_row,top_left_column\bottom_right_row,bottom_right_column| - already made
	// AB=askla{i}_OPT=Z2*AA2*$W$26*G2
	AB2=expxoptmx*Z[|.,.\rows(Z)-1,.|]:*AA[|.,.\rows(AA)-1,.|]:*G
	AB=(AB1\AB2)
	riskscoreopt=colsum(AB)
	riskvectoropt[i]=riskscoreopt

	//--------------------- Normal risk score ---------------------

	// Q - make input vector into column vector
	Xnor=X[1]\X[2]\ln(125)\0\0\ln(22.5)\0

	// T26 - Xnor*beta
	xnorbeta=Xnor'*coeffx
	dxnorbeta=Xnor'*coeffdx

	// $W$39=EXP(T48-W2)
	expxnormx=exp(xnorbeta-mxbeta)

	// $X$39=EXP(U48-X2)
	expdxnormdx=exp(dxnorbeta-mdxbeta)

	// building vector AB - askla{1}_OPT=$W$26*(-LN(E2)) and askla{i}_OPT=Z2*AA2*$W$26*G2
	// building vector AF - askla{1}_NOR=$W$39*(-LN(E2)) and askla{i}_NOR=AD2*AE2*$W$39*G2
	AF1=expxnormx*(-ln(E[1,1]))
	// AD - acs**exp(x-mx) NOR=POWER(E2,$W$39)
	AD=E:^expxnormx
	// AE - ads**exp(dx-mdx) NOR=POWER(J2,$X$39)
	AE=J:^expdxnormdx
	// G - log{i-1}-log{i}, subset |top_left_row,top_left_column\bottom_right_row,bottom_right_column| - already made
	//askla{i}_NOR=AD2*AE2*$W$39*G2
	AF2=expxnormx*AD[|.,.\rows(AD)-1,.|]:*AE[|.,.\rows(AE)-1,.|]:*G
	AF=(AF1\AF2)
	riskscorenor=colsum(AF)
	riskvectornor[i]=riskscorenor
	}

// export risk score vectors into pre-created Stata variables

st_store(., "hriskcvd30", touse, riskvector)
st_store(., "hriskcvd30opt", touse, riskvectoropt)
st_store(., "hriskcvd30nor", touse, riskvectornor)

}

end



//------------------ Mata sub-program computing full CVD risk (lipids) ------------------------------

mata

real matrix RISKLF (string scalar male, string scalar age, string scalar sbp, string scalar tcl, string scalar hdl, string scalar smoke, string scalar trtbp, string scalar diab, string scalar touse)
{
real scalar N, xbeta
real vector _sex, _logage, _logsbp, _logtcl, _loghdl, _smoke, _trtbp, _newdiab, _inputvars,  riskvector, riskvectoropt, riskvectornor

// create matrix which is a view on the current Stata dataset

st_view(inputvars=., ., (tokens(male), tokens(age), tokens(sbp), tokens(tcl), tokens(hdl), tokens(smoke), tokens(trtbp), tokens(diab)), touse)

// inputvars[1,.]

// Transform inputed values: sex, logage, logsbp, lntcl, lnhdl, smoke, trtbp, newdiab

_sex=inputvars[.,1]
_logage=ln(inputvars[.,2])
_logsbp=ln(inputvars[.,3])
_logtcl=ln(inputvars[.,4])
_loghdl=ln(inputvars[.,5])
_smoke=inputvars[.,6]
_trtbp=inputvars[.,7]
_newdiab=inputvars[.,8]

_inputvars=_sex,_logage,_logsbp,_logtcl,_loghdl,_smoke,_trtbp,_newdiab

//_inputvars[1,.]

// Load initial data (it's loaded in pieces due to Stata memory limitation)

// Load data for normal CVD risk

#delimit;

E1=(0.999874296\0.999874296\0.999874296\0.999748519\0.999622704\0.999622704\0.999496803\0.999370849\0.999370849\0.999370849\0.999244705\0.999244705\0.99911853\0.99911853\0.99911853\0.998992174\0.998865757\0.998739283\0.998612704\0.998612704\0.998612704\0.998486005\0.998359261\0.998232504\0.998105675\0.997978839\0.997851933\0.997851933\0.997851933\0.997724749\
0.997597515\0.997470181\0.997342824\0.997342824\0.997215427\0.997215427\0.997087879\0.996960256\0.996960256\0.996960256\0.996960256\0.996832045\0.996703795\0.996575514\0.996447164\0.996318743\0.996318743\0.996190148\0.996061122\0.995932052\0.995802901\0.995673626\0.995544206\0.995414666\0.995284937\0.99515505\0.995025081\0.995025081\0.995025081\0.994894865\
0.994764513\0.994634039\0.994503411\0.994503411\0.994503411\0.994372668\0.994241919\0.994110783\0.993979643\0.99384838\0.993717032\0.99358538\0.993453732\0.993322039\0.993190335\0.993058573\0.992926616\0.992794617\0.992662602\0.992530547\0.992398482\0.992266422\0.99213433\0.992002235\0.991870046\0.991737769\0.991605443\0.991473096\0.991340627\0.991208044\
0.991075452\0.990942808\0.990810007\0.990677101\0.990544031\0.99041084\0.990277491\0.990144063\0.990010525\0.989877001\0.989743334\0.989609467\0.989475554\0.989341599\0.989207647\0.989073673\0.988939555\0.988805391\0.988671231\0.988671231\0.988537062\0.988402883\0.988268691\0.988134423\0.988000123\0.987865777\0.987731408\0.987596881\0.98746234\0.987327798\
0.987193229\0.987058645\0.986924022\0.986789375\0.986654394\0.986519398\0.986384286\0.986384286\0.986249068\0.986113843\0.986113843\0.985978289\0.985978289\0.985978289\0.985842542\0.985706568\0.985570488\0.985434336\0.985297906\0.985297906\0.985161383\0.985024827\0.984888162\0.984751508\0.984614825\0.98447809\0.984341371\0.984341371\0.984204554\0.984067557\
0.983930384\0.983930384\0.983793063\0.983655684\0.983518256\0.983380646\0.983242918\0.983105116\0.982967179\0.982967179\0.982967179\0.982967179\0.982967179\0.982828783\0.982690349\0.982551894\0.982412785\0.982273637\0.982134487\0.981995136\0.981855744\0.981716262\0.981576767\0.981576767\0.981576767\0.981576767\0.981437115\0.981437115\0.981297367\0.981157515\
0.981017606\0.980877662\0.980737695\0.980597721\0.980457516\0.980317278\0.980317278\0.980176964\0.98003653\0.979896053\0.979755489\0.979614885\0.979474273\0.979333587\0.979192752\0.979051888\0.978910968\0.978769968\0.978628829\0.978628829\0.978628829\0.978487594\0.978346298\0.978346298\0.978204743\0.978204743\0.978062894\0.978062894\0.977920908\0.977778907\
0.977778907\0.977636652\0.977494371\0.977352026\0.977352026\0.977209622\0.977066907\0.976924135\0.976924135\0.976781251\0.97663828\0.976495262\0.976352089\0.976208855\0.976065627\0.976065627\0.975922003\0.975778365\0.975634644\0.975490766\0.975346853\0.975202466\0.975202466\0.975057997\0.974913545\0.974769045\0.974624521\0.974624521\0.974479884\0.974335227);
E2=(0.974335227\0.974190549\0.974190549\0.974190549\0.974045614\0.973900631\0.973755587\0.973610487\0.973465252\0.973465252\0.973465252\0.973319553\0.973173871\0.973028141\0.9728824\0.9728824\0.972736535\0.972736535\0.972590658\0.972444751\0.972298724\0.972152695\0.972006611\0.971860537\0.971714419\0.971568153\0.971421695\0.971275181\0.971128648\0.970982126\
0.970835606\0.970835606\0.97068888\0.970542117\0.970395279\0.970395279\0.970395279\0.970248233\0.970101034\0.969953813\0.969806594\0.969659159\0.969511667\0.969511667\0.969364022\0.969364022\0.969216258\0.969068343\0.968920351\0.968772291\0.968772291\0.968624043\0.968475778\0.968327498\0.968179149\0.968030784\0.967882356\0.967882356\0.967733771\0.96758498\
0.967436187\0.967287403\0.967287403\0.967138554\0.966989689\0.966840812\0.966691942\0.966543078\0.966394154\0.966245095\0.966095981\0.965946816\0.965797602\0.965648227\0.965498855\0.965349408\0.965349408\0.965349408\0.965199718\0.965050012\0.965050012\0.964900264\0.964900264\0.964750366\0.964600352\0.964450344\0.964300261\0.964300261\0.964149989\0.963999688\
0.963849367\0.963698908\0.963548398\0.963548398\0.96339782\0.96339782\0.963247188\0.963096443\0.963096443\0.962945659\0.962945659\0.962945659\0.962794662\0.962643638\0.962492525\0.962492525\0.962341371\0.962190203\0.962039023\0.961887771\0.961736518\0.961585269\0.961585269\0.961433927\0.961433927\0.961282421\0.961130891\0.960979285\0.960979285\0.960827614\
0.960675856\0.960524041\0.960524041\0.960524041\0.960372039\0.960219955\0.960219955\0.96006774\0.959915468\0.959763116\0.959610664\0.959458171\0.959458171\0.959305516\0.959305516\0.959152768\0.959152768\0.958999926\0.958847066\0.958694091\0.958541094\0.958388101\0.958388101\0.958234856\0.958234856\0.958081437\0.958081437\0.957927899\0.957774117\0.957774117\
0.957774117\0.95762008\0.957466054\0.957466054\0.957311767\0.957157467\0.957003047\0.956848629\0.95669398\0.956539186\0.956539186\0.956539186\0.956384152\0.956384152\0.956228946\0.956073689\0.955918444\0.955763108\0.955607463\0.955607463\0.955451746\0.95529599\0.95529599\0.95529599\0.955139997\0.954983898\0.954983898\0.954983898\0.954827682\0.95467124\
0.95467124\0.95467124\0.954514566\0.954357871\0.954201081\0.954044289\0.953887436\0.953730557\0.953573636\0.953573636\0.953416489\0.953416489\0.953259332\0.953102147\0.95294493\0.95294493\0.952787098\0.952787098\0.952787098\0.95262901\0.952470904\0.952470904\0.952312712\0.952154381\0.952154381\0.951995979\0.951837124\0.951678243\0.951519262\0.951360094\
0.951200825\0.951041322\0.951041322\0.950881796\0.950722283\0.950562708\0.950403016\0.950243224\0.950083389\0.949923468\0.949763514\0.949603492\0.949443324\0.949283138\0.949283138\0.949122859\0.948962447\0.948962447\0.948962447\0.948801463\0.948801463\0.948801463\0.948801463\0.948640018\0.9484783\0.948316452\0.948154508\0.947992378\0.947830268\0.947668157);
E3=(0.947506027\0.947343857\0.947181604\0.947181604\0.947181604\0.947019218\0.946856719\0.946694133\0.946531507\0.946368822\0.946206111\0.946206111\0.946043173\0.945880246\0.945880246\0.945717253\0.945554279\0.945391222\0.945391222\0.945391222\0.94522795\0.94522795\0.94506463\0.94506463\0.944901254\0.944737448\0.944737448\0.944737448\0.944573344\0.944409121\
0.944244735\0.944080355\0.944080355\0.943915684\0.943750882\0.943586066\0.943421226\0.94325639\0.94325639\0.943091276\0.94292615\0.94276104\0.94276104\0.942595815\0.942430592\0.942265238\0.942265238\0.942265238\0.942265238\0.942098959\0.941932557\0.941765837\0.941599089\0.941432347\0.941432347\0.941265422\0.94109848\0.940931478\0.94076443\0.94059739\
0.940430182\0.940430182\0.94026287\0.940095474\0.940095474\0.939927979\0.939760479\0.939592951\0.939592951\0.939592951\0.939425277\0.939257533\0.939089674\0.939089674\0.938921714\0.938753565\0.938753565\0.938585394\0.938585394\0.938417032\0.938248634\0.938248634\0.938079961\0.938079961\0.938079961\0.938079961\0.937911071\0.937742136\0.937742136\0.93757305\
0.937403973\0.937403973\0.937234802\0.937234802\0.937065446\0.937065446\0.936895925\0.936895925\0.936726281\0.936556554\0.936556554\0.936556554\0.936556554\0.936556554\0.936386333\0.936386333\0.936215834\0.936215834\0.936045215\0.936045215\0.936045215\0.935874349\0.935703428\0.935532411\0.935361239\0.935361239\0.93518953\0.935017742\0.934845766\0.9346738\
0.934501798\0.934501798\0.934329615\0.934157337\0.933984956\0.93381257\0.93381257\0.933640143\0.933467678\0.93329518\0.933122656\0.932950062\0.932777266\0.932604444\0.932431572\0.93225867\0.932085681\0.931912703\0.931739719\0.931739719\0.931566341\0.931392848\0.931219264\0.931045673\0.931045673\0.930872022\0.930698337\0.930524646\0.930350819\0.930176954\
0.930003027\0.930003027\0.930003027\0.930003027\0.929828639\0.929828639\0.929828639\0.929828639\0.929652749\0.929652749\0.929476694\0.929476694\0.929476694\0.929300451\0.929124189\0.928947889\0.928771582\0.928595204\0.928418663\0.928242028\0.928242028\0.928065115\0.927888069\0.927888069\0.92771081\0.927533551\0.927356318\0.927179096\0.927001895\0.926824608\
0.926646967\0.926469337\0.926469337\0.92629152\0.926113674\0.925935725\0.925757649\0.925757649\0.925578986\0.925400302\0.925221537\0.925042706\0.924863552\0.924684358\0.924684358\0.924504938\0.924325533\0.924146117\0.923966608\0.923787088\0.923607551\0.923427979\0.9232484\0.9232484\0.92306866\0.922888846\0.922888846\0.922708672\0.922528268\0.922347876\
0.922167375\0.922167375\0.921986856\0.921806159\0.921806159\0.921806159\0.92162507\0.92162507\0.92162507\0.92162507\0.92162507\0.921443234\0.921261361\0.921079197\0.921079197\0.920896783\0.920896783\0.920896783\0.920896783\0.920713653\0.920713653\0.920713653\0.920530377\0.920346958\0.920346958\0.920346958\0.920346958\0.920163166\0.920163166\0.919979081);
E4=(0.919795018\0.919795018\0.919609774\0.919609774\0.919424428\0.919239076\0.919053671\0.918868201\0.918682752\0.918682752\0.918682752\0.918682752\0.918682752\0.918497039\0.918497039\0.918497039\0.918497039\0.918497039\0.9183108\0.91812455\0.91812455\0.91812455\0.917937621\0.917750694\0.917563749\0.917376753\0.917189711\0.917189711\0.917002559\0.917002559\
0.916815366\0.916815366\0.916815366\0.916627824\0.916627824\0.916440088\0.916252234\0.916064351\0.916064351\0.915876269\0.915876269\0.915688078\0.915688078\0.915499827\0.915311568\0.915123275\0.915123275\0.915123275\0.914934867\0.914746463\0.914558022\0.914369586\0.914369586\0.914180977\0.913992315\0.913803603\0.913614851\0.913426069\0.913426069\0.913237126\
0.913237126\0.913047569\0.912857865\0.912857865\0.91266803\0.91266803\0.91266803\0.912477752\0.912287388\0.912096627\0.911905882\0.911714997\0.911714997\0.911714997\0.911714997\0.911524024\0.911332997\0.911141857\0.911141857\0.910950685\0.910759508\0.910759508\0.910568264\0.910568264\0.910568264\0.910568264\0.910376674\0.910184975\0.909993113\0.909993113\
0.909800823\0.909608489\0.909416163\0.909223778\0.909031271\0.908838765\0.908646183\0.908453521\0.908453521\0.908453521\0.908260642\0.908260642\0.908067664\0.907874683\0.907681695\0.907488697\0.907295651\0.907102553\0.906909453\0.906716241\0.906523016\0.906523016\0.90632935\0.90632935\0.906135437\0.906135437\0.90594136\0.90594136\0.905746952\0.905552539\
0.905358072\0.905358072\0.905358072\0.905163472\0.904968846\0.904774128\0.904579394\0.904579394\0.9043846\0.904189751\0.903994798\0.903799878\0.903604956\0.903409939\0.903409939\0.903409939\0.903409939\0.903409939\0.903213881\0.903017815\0.903017815\0.903017815\0.902820915\0.902820915\0.902624009\0.902427131\0.902230219\0.902230219\0.902230219\0.902033088\
0.901835698\0.901638213\0.901440706\0.901243196\0.901045538\0.900847777\0.900649884\0.900452019\0.900452019\0.900254121\0.900254121\0.900055613\0.900055613\0.899857064\0.89965846\0.89965846\0.89965846\0.89965846\0.899459566\0.899260645\0.899260645\0.899061532\0.898862287\0.898662896\0.898463488\0.898264076\0.898264076\0.898264076\0.898064172\0.897864135\
0.897663991\0.897463781\0.897263576\0.897263576\0.897063234\0.896862904\0.896862904\0.896862904\0.896662148\0.896662148\0.896461137\0.896260064\0.896058867\0.895857675\0.895857675\0.895656406\0.895656406\0.895455126\0.895253674\0.895052019\0.894850224\0.894850224\0.894648221\0.894446117\0.894243613\0.894243613\0.894040779\0.893837918\0.893634982\0.893431979\
0.89322866\0.893025341\0.893025341\0.893025341\0.893025341\0.892821584\0.892821584\0.892617657\0.89241359\0.892209494\0.892005322\0.891801121\0.891596892\0.891596892\0.891392109\0.891392109\0.891187281\0.890982391\0.890982391\0.890777347\0.890572209\0.890572209\0.890366852\0.890161347\0.890161347\0.889955602\0.889955602\0.889749803\0.889749803\0.889543939);
E5=(0.889543939\0.889337758\0.889337758\0.889131262\0.888924724\0.888924724\0.888717971\0.888511135\0.888511135\0.888304099\0.888097049\0.887890021\0.887683001\0.887475902\0.887268644\0.887268644\0.887060923\0.886852859\0.886644717\0.886644717\0.886436368\0.886227689\0.886018586\0.885809417\0.885600209\0.885390993\0.885390993\0.885390993\0.885390993\0.885390993\
0.885390993\0.885181294\0.884971544\0.884971544\0.884761771\0.884551852\0.884341684\0.884131463\0.884131463\0.884131463\0.883920966\0.883920966\0.883710289\0.883499612\0.883288923\0.883288923\0.883288923\0.883077797\0.882866575\0.882655217\0.882655217\0.882655217\0.882443692\0.882443692\0.882232127\0.882020582\0.881809046\0.881809046\0.881597278\0.881597278\
0.881597278\0.881597278\0.881385332\0.881385332\0.881173244\0.881173244\0.880960961\0.880748449\0.880535932\0.880323455\0.880323455\0.880110897\0.879898298\0.879685403\0.879472504\0.879259551\0.879046597\0.879046597\0.878833583\0.878620435\0.878620435\0.878407089\0.878193528\0.87797989\0.87797989\0.877766106\0.877766106\0.877552255\0.877338394\0.877124474\
0.876910584\0.876910584\0.876696386\0.876482236\0.876482236\0.876267949\0.876267949\0.876267949\0.876052993\0.876052993\0.876052993\0.875837798\0.875622469\0.875406998\0.87519146\0.874975701\0.8747599\0.8747599\0.874543657\0.874543657\0.874543657\0.874327144\0.874327144\0.874327144\0.874110084\0.873892955\0.873675807\0.873458517\0.873458517\0.873240987\
0.873240987\0.873023302\0.872805522\0.872587648\0.872587648\0.872587648\0.872587648\0.872368834\0.872150049\0.872150049\0.872150049\0.871931106\0.871712197\0.871712197\0.871493218\0.871274064\0.871054596\0.871054596\0.8708348\0.8708348\0.870614949\0.870614949\0.870394773\0.870174574\0.869954291\0.869954291\0.86973362\0.869512979\0.869512979\0.86929231\
0.869071343\0.868850261\0.86862915\0.868407512\0.868185836\0.867964037\0.867742082\0.867520172\0.867298213\0.867298213\0.867298213\0.867075628\0.866852697\0.866852697\0.866629677\0.866406566\0.866183424\0.866183424\0.865960023\0.865960023\0.865960023\0.865736148\0.865736148\0.865736148\0.865512062\0.865512062\0.865512062\0.865287049\0.865062002\0.864836735\
0.864611487\0.864611487\0.864385729\0.864385729\0.864385729\0.86415976\0.863933685\0.863707555\0.86348119\0.863254855\0.863254855\0.863254855\0.863028284\0.862801682\0.862575048\0.862575048\0.862347989\0.86212096\0.861893882\0.861666802\0.861439573\0.861212112\0.861212112\0.861212112\0.861212112\0.861212112\0.861212112\0.861212112\0.860982893\0.860753706\
0.860524541\0.860295372\0.860295372\0.860065978\0.860065978\0.859836596\0.859607114\0.859377453\0.859147699\0.859147699\0.859147699\0.858917365\0.858687072\0.858456798\0.858226511\0.858226511\0.858226511\0.858226511\0.858226511\0.857995176\0.857995176\0.857995176\0.857995176\0.857763435\0.857763435\0.857531462\0.857531462\0.857299304\0.857299304\0.857299304);
E6=(0.857066579\0.85683364\0.85683364\0.856600282\0.856600282\0.856600282\0.856366617\0.856366617\0.856366617\0.856132586\0.856132586\0.855898226\0.85566377\0.855429129\0.855194455\0.855194455\0.854959606\0.854959606\0.854724719\0.854489206\0.854253516\0.854017711\0.853781849\0.853545799\0.853545799\0.853309622\0.853073431\0.852837286\0.852600667\0.852363798\
0.852363798\0.852363798\0.852126598\0.852126598\0.851889371\0.851889371\0.851652044\0.851414353\0.851414353\0.851414353\0.85117647\0.850938153\0.850699695\0.850460943\0.850222188\0.850222188\0.850222188\0.850222188\0.850222188\0.849982036\0.849982036\0.849982036\0.849740747\0.849499468\0.849499468\0.849257957\0.849016442\0.848774496\0.848774496\0.848774496\
0.848531656\0.848288833\0.848045654\0.848045654\0.847802116\0.847802116\0.847802116\0.84755819\0.84731423\0.84731423\0.84731423\0.847070074\0.846825718\0.846825718\0.846581211\0.846336635\0.846336635\0.846091775\0.845846899\0.845601658\0.845356466\0.845356466\0.845110502\0.844864403\0.844618154\0.844371932\0.844125718\0.844125718\0.843879323\0.843879323\
0.843632735\0.843386015\0.84313927\0.842892503\0.842645168\0.842397757\0.842397757\0.842397757\0.842149869\0.842149869\0.842149869\0.841901255\0.841901255\0.841901255\0.841651628\0.841651628\0.8414013\0.841150794\0.840899408\0.840647476\0.840395289\0.840142876\0.839890232\0.839890232\0.839890232\0.839890232\0.83963548\0.83963548\0.839380292\0.839380292\
0.839124777\0.838869302\0.838869302\0.838869302\0.838869302\0.838612298\0.838354766\0.838354766\0.838354766\0.838354766\0.838095497\0.838095497\0.837835339\0.837835339\0.837574313\0.837574313\0.837312793\0.837312793\0.837312793\0.837312793);

J1=(1\0.999874446\0.999748865\0.999748865\0.999748865\0.999623134\0.999623134\0.999623134\0.999497251\0.999371268\0.999371268\0.9992452\0.9992452\0.999119087\0.998992881\0.998992881\0.998992881\0.998992881\0.998992881\0.998866203\0.998739535\0.998739535\0.998739535\0.998739535\0.998739535\0.998739535\0.998739535\0.998612524\0.998485376\0.998485376\
0.998485376\0.998485376\0.998485376\0.998357988\0.998357988\0.998230522\0.998230522\0.998230522\0.99810285\0.997974895\0.997846947\0.997846947\0.997846947\0.997846947\0.997846947\0.997846947\0.997718564\0.997718564\0.997718564\0.997718564\0.997718564\0.997718564\0.997718564\0.997718564\0.997718564\0.997718564\0.997718564\0.997588854\0.997459127\0.997459127\
0.997459127\0.997459127\0.997459127\0.997328722\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\
0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997198313\0.997064238\0.997064238\0.997064238\0.997064238\0.997064238\0.997064238\0.997064238\0.997064238\0.997064238\0.997064238\0.997064238\
0.997064238\0.997064238\0.997064238\0.997064238\0.997064238\0.997064238\0.997064238\0.996929035\0.996929035\0.996929035\0.996793696\0.996793696\0.996658004\0.996522198\0.996522198\0.996522198\0.996522198\0.996522198\0.996522198\0.996385463\0.996385463\0.996385463\0.996385463\0.996385463\0.996385463\0.996385463\0.996385463\0.996248477\0.996248477\0.996248477\
0.996248477\0.996111222\0.996111222\0.996111222\0.996111222\0.996111222\0.996111222\0.996111222\0.996111222\0.99597303\0.995834834\0.995696584\0.995558221\0.995558221\0.995558221\0.995558221\0.995558221\0.995558221\0.995558221\0.995558221\0.995558221\0.995558221\0.995558221\0.995418982\0.995279727\0.995140439\0.995140439\0.995001025\0.995001025\0.995001025\
0.995001025\0.995001025\0.995001025\0.995001025\0.995001025\0.995001025\0.994861076\0.994861076\0.994861076\0.994861076\0.994861076\0.994861076\0.994861076\0.994861076\0.994861076\0.994861076\0.994861076\0.994861076\0.994861076\0.994720088\0.994579115\0.994579115\0.994579115\0.994437815\0.994437815\0.994296454\0.994296454\0.9941548\0.9941548\0.9941548\
0.994012805\0.994012805\0.994012805\0.994012805\0.993870559\0.993870559\0.993870559\0.993870559\0.993727889\0.993727889\0.993727889\0.993727889\0.993727889\0.993727889\0.993727889\0.99358467\0.99358467\0.99358467\0.99358467\0.99358467\0.99358467\0.99358467\0.993440487\0.993440487\0.993440487\0.993440487\0.993440487\0.993296187\0.993296187\0.99329618);
J2=(0.993151759\0.993151759\0.993007258\0.992862583\0.992862583\0.992862583\0.992862583\0.992862583\0.992862583\0.992717442\0.992571828\0.992571828\0.992571828\0.992571828\0.992571828\0.992425943\0.992425943\0.992279943\0.992279943\0.992279943\0.992279943\0.992279943\0.992279943\0.992279943\0.992279943\0.992279943\0.992279943\0.992279943\0.992279943\0.992279943\
0.992279943\0.992133146\0.992133146\0.992133146\0.992133146\0.991986093\0.991838991\0.991838991\0.991838991\0.991838991\0.991838991\0.991838991\0.991838991\0.991691343\0.991691343\0.991543489\0.991543489\0.991543489\0.991543489\0.991543489\0.991395157\0.991395157\0.991395157\0.991395157\0.991395157\0.991395157\0.991395157\0.99124619\0.99124619\0.99124619\
0.99124619\0.99124619\0.991096646\0.991096646\0.991096646\0.991096646\0.991096646\0.991096646\0.991096646\0.991096646\0.991096646\0.991096646\0.991096646\0.991096646\0.991096646\0.991096646\0.990946141\0.990795565\0.990795565\0.990795565\0.990644811\0.990644811\0.990494032\0.990494032\0.990494032\0.990494032\0.990494032\0.990342859\0.990342859\0.990342859\
0.990342859\0.990342859\0.990342859\0.990191349\0.990191349\0.990039739\0.990039739\0.990039739\0.989887945\0.989887945\0.989736103\0.989584273\0.989584273\0.989584273\0.989584273\0.989432141\0.989432141\0.989432141\0.989432141\0.989432141\0.989432141\0.989432141\0.989279712\0.989279712\0.989127204\0.989127204\0.989127204\0.989127204\0.988974423\0.988974423\
0.988974423\0.988974423\0.988821336\0.988668143\0.988668143\0.988668143\0.988514805\0.988514805\0.988514805\0.988514805\0.988514805\0.988514805\0.988360876\0.988360876\0.988206737\0.988206737\0.988052415\0.988052415\0.988052415\0.988052415\0.988052415\0.988052415\0.987897673\0.987897673\0.987742785\0.987742785\0.987587711\0.987587711\0.987587711\0.987431977\
0.987276179\0.987276179\0.987276179\0.987120078\0.987120078\0.987120078\0.987120078\0.987120078\0.987120078\0.987120078\0.986963378\0.986806671\0.986806671\0.986649751\0.986649751\0.986649751\0.986649751\0.986649751\0.986649751\0.986491864\0.986491864\0.986491864\0.986333708\0.986175447\0.986175447\0.986175447\0.986016981\0.985858436\0.985858436\0.985858436\
0.985699689\0.985540838\0.985540838\0.985540838\0.985540838\0.985540838\0.985540838\0.985540838\0.985540838\0.98538154\0.98538154\0.985222124\0.985222124\0.985222124\0.985222124\0.985062358\0.985062358\0.984902366\0.984742339\0.984742339\0.984742339\0.984582229\0.984582229\0.984582229\0.984421909\0.984421909\0.984421909\0.984421909\0.984421909\0.984421909\
0.984421909\0.984421909\0.984259986\0.984259986\0.984259986\0.984259986\0.984259986\0.984259986\0.984259986\0.984259986\0.984259986\0.984259986\0.984259986\0.984259986\0.98409695\0.98409695\0.98409695\0.983933535\0.98376996\0.98376996\0.983606107\0.983442194\0.983278228\0.983278228\0.983278228\0.983278228\0.983278228\0.983278228\0.983278228\0.98327822);
J3=(0.983278228\0.983278228\0.983278228\0.983113004\0.982947774\0.982947774\0.982947774\0.982947774\0.982947774\0.982947774\0.982947774\0.982781856\0.982781856\0.982781856\0.982615818\0.982615818\0.982615818\0.982615818\0.982449571\0.982283314\0.982283314\0.982116959\0.982116959\0.981950557\0.981950557\0.981950557\0.981783809\0.981616841\0.981616841\0.981616841\
0.981616841\0.981616841\0.981449099\0.981449099\0.981449099\0.981449099\0.981449099\0.981449099\0.981280862\0.981280862\0.981280862\0.981280862\0.98111221\0.98111221\0.98111221\0.98111221\0.980943124\0.980773928\0.980604674\0.980604674\0.980604674\0.980604674\0.980604674\0.980604674\0.980434768\0.980434768\0.980434768\0.980434768\0.980434768\0.980434768\
0.980434768\0.980264418\0.980264418\0.980264418\0.980093896\0.980093896\0.980093896\0.980093896\0.979923115\0.979752338\0.979752338\0.979752338\0.979752338\0.97958123\0.97958123\0.97958123\0.979409962\0.979409962\0.979238546\0.979238546\0.979238546\0.979066823\0.979066823\0.978894732\0.978722655\0.978550582\0.978550582\0.978550582\0.978378325\0.978378325\
0.978378325\0.978206037\0.978206037\0.978033489\0.978033489\0.977860739\0.977860739\0.977687929\0.977687929\0.977687929\0.977514828\0.977341681\0.977168541\0.976995316\0.976995316\0.976821653\0.976821653\0.976647876\0.976647876\0.976473881\0.976299881\0.976299881\0.976299881\0.976299881\0.976299881\0.97612507\0.97612507\0.97612507\0.97612507\0.97612507\
0.97612507\0.975949811\0.975949811\0.975949811\0.975949811\0.975949811\0.975774192\0.975774192\0.975774192\0.975774192\0.975774192\0.975774192\0.975774192\0.975774192\0.975774192\0.975774192\0.975774192\0.975774192\0.975774192\0.975597025\0.975597025\0.975597025\0.975597025\0.975597025\0.975419464\0.975419464\0.975419464\0.975419464\0.975419464\0.975419464\
0.975419464\0.975241279\0.97506285\0.974884412\0.974884412\0.974705183\0.974525779\0.974346036\0.974346036\0.974166058\0.974166058\0.973985978\0.973805892\0.973805892\0.973805892\0.973805892\0.973805892\0.973805892\0.973805892\0.973805892\0.973624804\0.973624804\0.973624804\0.973443533\0.973443533\0.973443533\0.973443533\0.973443533\0.973443533\0.973443533\
0.973443533\0.973443533\0.973261196\0.973261196\0.973261196\0.973261196\0.973261196\0.973077915\0.973077915\0.973077915\0.973077915\0.973077915\0.973077915\0.973077915\0.972893467\0.972893467\0.972893467\0.972893467\0.972893467\0.972893467\0.972893467\0.972893467\0.972893467\0.972708531\0.972708531\0.972708531\0.972523005\0.972523005\0.972523005\0.972523005\
0.972523005\0.972337159\0.972337159\0.972337159\0.972151053\0.971964795\0.971964795\0.971778277\0.97159176\0.971405055\0.97121817\0.97121817\0.97121817\0.97121817\0.971030539\0.971030539\0.970842602\0.970654352\0.970466015\0.970466015\0.970277437\0.970088742\0.970088742\0.970088742\0.969899877\0.969710995\0.969522042\0.969522042\0.969332687\0.96933268);
J4=(0.969332687\0.96914297\0.96914297\0.968952631\0.968952631\0.968952631\0.968952631\0.968952631\0.968952631\0.968761902\0.968571196\0.968380452\0.968189699\0.968189699\0.96799888\0.967808063\0.967617117\0.96742616\0.96742616\0.96742616\0.967234784\0.967043239\0.967043239\0.967043239\0.967043239\0.967043239\0.967043239\0.966851234\0.966851234\0.966659048\
0.966659048\0.966466853\0.966274658\0.966274658\0.96608215\0.96608215\0.96608215\0.96608215\0.965889231\0.965889231\0.965696071\0.965696071\0.965502857\0.965502857\0.965502857\0.965502857\0.965309494\0.965116135\0.965116135\0.965116135\0.965116135\0.965116135\0.964922428\0.964922428\0.964922428\0.964922428\0.964922428\0.964922428\0.964728309\0.964728309\
0.964534053\0.964534053\0.964534053\0.964339069\0.964339069\0.964143914\0.963948478\0.963948478\0.963948478\0.963948478\0.963948478\0.963948478\0.963751942\0.963555389\0.963358867\0.963358867\0.963358867\0.963358867\0.963162087\0.963162087\0.963162087\0.962965283\0.962965283\0.962768304\0.962571324\0.96237436\0.96237436\0.96237436\0.96237436\0.96217668\
0.96217668\0.96217668\0.96217668\0.96217668\0.96217668\0.96217668\0.96217668\0.96217668\0.961977764\0.961778794\0.961778794\0.961579764\0.961579764\0.961579764\0.961579764\0.961579764\0.961579764\0.961579764\0.961579764\0.961579764\0.961579764\0.961379915\0.961379915\0.961179561\0.961179561\0.960979108\0.960979108\0.96077856\0.96077856\0.96077856\
0.96077856\0.960577673\0.960376757\0.960376757\0.960376757\0.960376757\0.960376757\0.960175407\0.960175407\0.960175407\0.960175407\0.960175407\0.960175407\0.960175407\0.959973554\0.959771348\0.959569041\0.95936666\0.95936666\0.95936666\0.959163843\0.958960561\0.958960561\0.958757238\0.958757238\0.958757238\0.958757238\0.958553777\0.958350338\0.958350338\
0.958350338\0.958350338\0.958350338\0.958350338\0.958350338\0.958350338\0.958350338\0.958350338\0.95814574\0.95814574\0.95794059\0.95794059\0.957735439\0.957735439\0.957735439\0.957529917\0.957324365\0.957118848\0.957118848\0.957118848\0.956913216\0.956913216\0.956913216\0.956913216\0.956913216\0.956913216\0.956706761\0.956500284\0.956500284\0.956500284\
0.956500284\0.956500284\0.956500284\0.956293226\0.956293226\0.956293226\0.956086006\0.955878548\0.955878548\0.955670863\0.955670863\0.955670863\0.955670863\0.955670863\0.955462751\0.955462751\0.955254592\0.955254592\0.955254592\0.955254592\0.955254592\0.955045911\0.955045911\0.955045911\0.955045911\0.954836577\0.954836577\0.954836577\0.954836577\0.954836577\
0.954836577\0.954836577\0.954626638\0.95441644\0.954206227\0.954206227\0.953995819\0.953995819\0.953995819\0.953995819\0.953995819\0.953995819\0.953995819\0.953784317\0.953784317\0.953572696\0.953572696\0.953572696\0.953360938\0.953360938\0.953360938\0.95314887\0.95314887\0.95314887\0.952936483\0.952936483\0.952723837\0.952723837\0.952511053\0.95251105);
J5=(0.952298089\0.952298089\0.952084634\0.952084634\0.952084634\0.951870849\0.951870849\0.951870849\0.951656776\0.951656776\0.951656776\0.951656776\0.951656776\0.951656776\0.951656776\0.951441981\0.951441981\0.951441981\0.951441981\0.951226288\0.951226288\0.951226288\0.951226288\0.951226288\0.951226288\0.951226288\0.951009621\0.950792957\0.950576114\0.950359254\
0.950142322\0.950142322\0.950142322\0.949925312\0.949925312\0.949925312\0.949925312\0.949925312\0.949707415\0.94948938\0.94948938\0.949271356\0.949271356\0.949271356\0.949271356\0.949053037\0.948834718\0.948834718\0.948834718\0.948834718\0.948616045\0.948397409\0.948397409\0.948178589\0.948178589\0.948178589\0.948178589\0.947959579\0.947959579\0.947740406\
0.947521175\0.947301961\0.947301961\0.947082717\0.947082717\0.946863079\0.946863079\0.946863079\0.946863079\0.946863079\0.946642963\0.946642963\0.946642963\0.946642963\0.946642963\0.946642963\0.946642963\0.946422233\0.946422233\0.946422233\0.946201141\0.946201141\0.946201141\0.946201141\0.945979522\0.945979522\0.945757694\0.945757694\0.945757694\0.945757694\
0.945757694\0.945535184\0.945535184\0.945535184\0.945312586\0.945312586\0.945089327\0.94486593\0.94486593\0.944642371\0.944418737\0.944418737\0.944418737\0.944418737\0.944418737\0.944418737\0.944418737\0.944194267\0.944194267\0.943969634\0.943744828\0.943744828\0.943519825\0.94329444\0.94329444\0.94329444\0.94329444\0.94329444\0.943068292\0.943068292\
0.942841957\0.942841957\0.942841957\0.942841957\0.942615177\0.942387889\0.942160388\0.942160388\0.942160388\0.941932711\0.941705059\0.941705059\0.941705059\0.941477211\0.941477211\0.941477211\0.941477211\0.941248607\0.941248607\0.941019614\0.941019614\0.940790577\0.940790577\0.940790577\0.940790577\0.940560597\0.940560597\0.940560597\0.940330397\0.940330397\
0.940330397\0.940330397\0.940330397\0.940330397\0.940330397\0.940330397\0.940330397\0.940330397\0.940330397\0.940098494\0.939866316\0.939866316\0.939866316\0.93963297\0.93963297\0.93963297\0.93963297\0.939399298\0.939399298\0.939165569\0.938931807\0.938931807\0.938697558\0.938463329\0.938463329\0.938228524\0.937993214\0.937993214\0.937993214\0.937993214\
0.937993214\0.937757176\0.937757176\0.937521102\0.93728504\0.93728504\0.93728504\0.93728504\0.93728504\0.93728504\0.937048371\0.936811724\0.936811724\0.936811724\0.936811724\0.936574582\0.936574582\0.936574582\0.936574582\0.936574582\0.936574582\0.936574582\0.936336544\0.936098525\0.935860329\0.935621421\0.935382522\0.935143483\0.935143483\0.935143483\
0.935143483\0.935143483\0.934904165\0.934904165\0.934664771\0.934664771\0.934664771\0.934664771\0.934664771\0.934424598\0.934184058\0.934184058\0.934184058\0.934184058\0.934184058\0.933943239\0.933702435\0.933461504\0.933220601\0.933220601\0.932979085\0.932737521\0.932495844\0.932495844\0.932254008\0.932254008\0.932011803\0.932011803\0.931769433\0.93152699);
J6=(0.931526992\0.931526992\0.931284152\0.931284152\0.93104114\0.930797877\0.930797877\0.93055433\0.93031076\0.93031076\0.930066837\0.930066837\0.930066837\0.930066837\0.930066837\0.929821942\0.929821942\0.929576968\0.929576968\0.929576968\0.929576968\0.929576968\0.929576968\0.929576968\0.929330349\0.929330349\0.929330349\0.929330349\0.929330349\0.929330349\
0.929082907\0.92883534\0.92883534\0.928587653\0.928587653\0.928339838\0.928339838\0.928339838\0.928091251\0.927842361\0.927842361\0.927842361\0.927842361\0.927842361\0.927842361\0.927591611\0.927340526\0.927089402\0.926838164\0.926838164\0.926586073\0.926333311\0.926333311\0.926333311\0.926080172\0.926080172\0.926080172\0.926080172\0.925825973\0.92557159\
0.92557159\0.92557159\0.92557159\0.925316333\0.925316333\0.925060562\0.924804746\0.924804746\0.924804746\0.924548653\0.924292574\0.924292574\0.924292574\0.924036176\0.924036176\0.924036176\0.92377903\0.92377903\0.92377903\0.92377903\0.92377903\0.92352153\0.92352153\0.92352153\0.92352153\0.92352153\0.92352153\0.923263092\0.923263092\0.923004562\
0.923004562\0.923004562\0.923004562\0.923004562\0.923004562\0.923004562\0.922744784\0.922484685\0.922484685\0.922224233\0.92196328\0.92196328\0.921701842\0.921439956\0.921439956\0.921177434\0.921177434\0.921177434\0.921177434\0.921177434\0.921177434\0.921177434\0.921177434\0.920910625\0.920643427\0.920376008\0.920376008\0.920107433\0.920107433\0.91983871\
0.91983871\0.91983871\0.91956953\0.919300087\0.919029826\0.919029826\0.919029826\0.918758565\0.918485908\0.918213173\0.918213173\0.917939832\0.917939832\0.917665669\0.917665669\0.917391013\0.917391013\0.917115523\0.916840007\0.91656423);

#delimit cr

// assemble pieces of initial data into single vectors

E=E1\E2\E3\E4\E5\E6
J=J1\J2\J3\J4\J5\J6

// Create matrix G (used below) - log{i-1}-log{i}, subset |top_left_row,top_left_column\bottom_right_row,bottom_right_column|
E1=E[|1,.\rows(E)-1,.|]
E2=E[|2,.\rows(E),.|]
G=ln(E1)-ln(E2)

// input initial values of parameters

coeffx=(0.34362\2.63588\1.8803\1.12673\-0.90941\0.59397\0.5232\0.68602)
coeffdx=(0.48123\3.39222\1.39862\-0.00439\0.16081\0.99858\0.19035\0.49756)
mxbeta=21.2932661180837
mdxbeta=20.1284069843952

// Define new vector containing risk scores
riskvector=J(rows(_inputvars),1,.)
riskvectoropt=J(rows(_inputvars),1,.)
riskvectornor=J(rows(_inputvars),1,.)

// loop over observations and compute risk score for each one

for (i=1; i<=rows(_inputvars); i++) {

	// select observation i of input data
	Xrow=_inputvars[| i,1\i,.|]

	//--------------------- Risk score ---------------------

	// Q - make input vector into column vector
	X=Xrow'

	// R - X*beta
	xbeta=X'*coeffx
	dxbeta=X'*coeffdx

	// $W$5
	expxmx=exp(xbeta-mxbeta)

	// $X$5=EXP(U10-X2)
	expdxmdx=exp(dxbeta-mdxbeta)

	// building vector M - askla{i}==$W$5*(-LN(E2)) and askla{i}=F2*K2*$W$5*G2
	M1=expxmx*(-ln(E[1,1]))

	// F - acs**exp(x-mx)=POWER(E2,$W$5)
	F=E:^expxmx

	// K - ads**exp(dx-mdx)=POWER(J3,$X$5)
	K=J:^expdxmdx

	// M=askla{i}=F2*K2*$W$5*G2
	M2=expxmx*F[|.,.\rows(F)-1,.|]:*K[|.,.\rows(K)-1,.|]:*G
	M=(M1\M2)

	riskscore=colsum(M)
	riskvector[i]=riskscore

	//--------------------- Optimal risk score ---------------------

	// Q - make input vector into column vector
	Xopt=X[1]\X[2]\ln(110)\ln(160)\ln(60)\0\0\0

	// T26 - Xopt*beta
	xoptbeta=Xopt'*coeffx
	dxoptbeta=Xopt'*coeffdx

	// $W$26=EXP(T35-W2)
	expxoptmx=exp(xoptbeta-mxbeta)

	// $X$26=EXP(U35-X2)
	expdxoptmdx=exp(dxoptbeta-mdxbeta)

	// building vector AB - askla{1}_OPT=$W$26*(-LN(E2)) and askla{i}_OPT=Z2*AA2*$W$26*G2
	AB1=expxoptmx*(-ln(E[1,1]))
	// Z - acs**exp(x-mx) OPT=POWER(E2,$W$26)
	Z=E:^expxoptmx
	// AA - ads**exp(dx-mdx) OPT==POWER(J2,$X$26)
	AA=J:^expdxoptmdx
	// G - log{i-1}-log{i}, subset |top_left_row,top_left_column\bottom_right_row,bottom_right_column| - already made
	// AB=askla{i}_OPT=Z2*AA2*$W$26*G2
	AB2=expxoptmx*Z[|.,.\rows(Z)-1,.|]:*AA[|.,.\rows(AA)-1,.|]:*G
	AB=(AB1\AB2)
	riskscoreopt=colsum(AB)
	riskvectoropt[i]=riskscoreopt

	//--------------------- Normal risk score ---------------------

	// Q - make input vector into column vector
	Xnor=X[1]\X[2]\ln(125)\ln(180)\ln(45)\0\0\0

	// T26 - Xnor*beta
	xnorbeta=Xnor'*coeffx
	dxnorbeta=Xnor'*coeffdx

	// $W$39=EXP(T48-W2)
	expxnormx=exp(xnorbeta-mxbeta)

	// $X$39=EXP(U48-X2)
	expdxnormdx=exp(dxnorbeta-mdxbeta)

	// building vector AB - askla{1}_OPT=$W$26*(-LN(E2)) and askla{i}_OPT=Z2*AA2*$W$26*G2
	// building vector AF - askla{1}_NOR=$W$39*(-LN(E2)) and askla{i}_NOR=AD2*AE2*$W$39*G2
	AF1=expxnormx*(-ln(E[1,1]))
	// AD - acs**exp(x-mx) NOR=POWER(E2,$W$39)
	AD=E:^expxnormx
	// AE - ads**exp(dx-mdx) NOR=POWER(J2,$X$39)
	AE=J:^expdxnormdx
	// G - log{i-1}-log{i}, subset |top_left_row,top_left_column\bottom_right_row,bottom_right_column| - already made
	//askla{i}_NOR=AD2*AE2*$W$39*G2
	AF2=expxnormx*AD[|.,.\rows(AD)-1,.|]:*AE[|.,.\rows(AE)-1,.|]:*G
	AF=(AF1\AF2)
	riskscorenor=colsum(AF)
	riskvectornor[i]=riskscorenor
	}

// export risk score vectors into pre-created Stata variables

st_store(., "friskcvdl30", touse, riskvector)
st_store(., "friskcvdl30opt", touse, riskvectoropt)
st_store(., "friskcvdl30nor", touse, riskvectornor)

}

end

//---------------------- Mata sub-program computing hard CVD risk ------------------------------

mata

real matrix RISKLH (string scalar male, string scalar age, string scalar sbp, string scalar tcl, string scalar hdl, string scalar smoke, string scalar trtbp, string scalar diab, string scalar touse)
{
real scalar N, xbeta
real vector _sex, _logage, _logsbp, _logtcl, _loghdl, _smoke, _trtbp, _newdiab, _inputvars,  riskvector, riskvectoropt, riskvectornor

// create matrix which is a view on the current Stata dataset

st_view(inputvars=., ., (tokens(male), tokens(age), tokens(sbp), tokens(tcl), tokens(hdl), tokens(smoke), tokens(trtbp), tokens(diab)), touse)

// inputvars[1,.]

// Transform inputed values: sex, logage, logsbp, lntcl, lnhdl, smoke, trtbp, newdiab

_sex=inputvars[.,1]
_logage=ln(inputvars[.,2])
_logsbp=ln(inputvars[.,3])
_logtcl=ln(inputvars[.,4])
_loghdl=ln(inputvars[.,5])
_smoke=inputvars[.,6]
_trtbp=inputvars[.,7]
_newdiab=inputvars[.,8]

_inputvars=_sex,_logage,_logsbp,_logtcl,_loghdl,_smoke,_trtbp,_newdiab

//_inputvars[1,.]

// Load initial data (it's loaded in pieces due to Stata memory limitation)

// Load data for hard CVD risk

#delimit;

E1=(0.99989223\0.99989223\0.99989223\0.9997844\0.999676537\0.999676537\0.999568591\0.999568591\0.999568591\0.999460474\0.999460474\0.999352333\0.999352333\0.999352333\0.999244022\0.999135647\0.999027221\0.998918678\0.998918678\0.998918678\0.998810045\0.998810045\0.998701294\0.998592538\0.998592538\0.998592538\0.99848344\0.998374296\0.998265086\0.998265086\
0.998155839\0.998155839\0.99804647\0.99804647\0.99804647\0.99804647\0.997936874\0.997827233\0.997717558\0.997717558\0.997607689\0.997497246\0.997386721\0.997276046\0.997165328\0.997054488\0.996943489\0.996832339\0.996832339\0.996832339\0.996720949\0.996609432\0.996497798\0.996497798\0.996497798\0.996386054\0.996274187\0.996162208\0.996050183\0.995938163\
0.995826068\0.995713817\0.995601373\0.995488809\0.995376109\0.995263343\0.995150588\0.995037182\0.995037182\0.99492377\0.994810334\0.994696853\0.994583373\0.994583373\0.994469691\0.994469691\0.994355653\0.994355653\0.994355653\0.994355653\0.994241355\0.99412703\0.99401268\0.993898284\0.993898284\0.993898284\0.993783683\0.993668887\0.99355398\0.993439001\
0.993324005\0.993324005\0.993324005\0.993324005\0.993324005\0.993208547\0.993092982\0.992977392\0.992861801\0.992746193\0.992630573\0.992630573\0.992514877\0.992514877\0.992514877\0.992398562\0.992398562\0.992282173\0.992165739\0.992049285\0.992049285\0.99193276\0.991816122\0.991699452\0.991582702\0.991465914\0.991349122\0.991232255\0.991115362\0.99099844\
0.990881387\0.990881387\0.990764257\0.990764257\0.990646852\0.990529399\0.99041194\0.99041194\0.990294124\0.990176296\0.990176296\0.990058191\0.989939936\0.989939936\0.989821602\0.989703254\0.989703254\0.989584663\0.989466015\0.989466015\0.989347321\0.98922858\0.98922858\0.989109735\0.988990804\0.988871843\0.988752729\0.988633557\0.98851439\0.98851439\
0.988394742\0.988274939\0.988274939\0.988155063\0.988035168\0.988035168\0.987915174\0.987795163\0.987675055\0.987675055\0.987554926\0.987554926\0.987554926\0.987434591\0.98731421\0.987193768\0.987193768\0.987193768\0.987072817\0.986951684\0.986951684\0.986830404\0.986830404\0.986709003\0.986587546\0.986466095\0.986344481\0.986344481\0.986222715\0.986222715\
0.986222715\0.986100683\0.985978497\0.985978497\0.985856165\0.985856165\0.985733729\0.985611141\0.985488483\0.985488483\0.985365602\0.985242703\0.985119406\0.985119406\0.985119406\0.984995886\0.984872351\0.984872351\0.984872351\0.984748582\0.984748582\0.984748582\0.984624675\0.984624675\0.984500637\0.984376575\0.984252493\0.984128395\0.984004254\0.984004254\
0.983879905\0.983879905\0.983879905\0.983755502\0.983755502\0.983755502\0.983630861\0.983630861\0.98350618\0.983381484\0.983256779\0.983132\0.983007145\0.98288229\0.98288229\0.982757355\0.982632406\0.982507442\0.982507442\0.982382346\0.982382346\0.982257185\0.982131956\0.982131956\0.982131956\0.98200657\0.98200657\0.981881085\0.981755544\0.98162994);
E2=(0.981504259\0.981378543\0.981378543\0.981252697\0.981126841\0.981126841\0.981000856\0.981000856\0.981000856\0.980874447\0.980874447\0.980747819\0.980747819\0.980621066\0.980621066\0.980494204\0.980494204\0.980494204\0.980366965\0.980366965\0.980239486\0.980111906\0.979984325\0.979984325\0.979984325\0.979856453\0.979856453\0.979728434\0.979600368\0.979472232\
0.979472232\0.979344037\0.979215808\0.979215808\0.979215808\0.979215808\0.97908735\0.97908735\0.97908735\0.97908735\0.978958646\0.978829869\0.978701086\0.978572207\0.978443194\0.978443194\0.978314057\0.978314057\0.978184737\0.978184737\0.97805528\0.977925797\0.977925797\0.977795841\0.977795841\0.977795841\0.977665752\0.977665752\0.977665752\0.977535425\
0.977405024\0.977405024\0.977405024\0.97727435\0.977143685\0.977012951\0.976882155\0.976751179\0.976751179\0.976751179\0.976620114\0.976488931\0.976488931\0.976357314\0.976225666\0.976225666\0.976093954\0.976093954\0.976093954\0.976093954\0.975961824\0.975829493\0.975696823\0.975564136\0.975431418\0.975431418\0.975431418\0.975298559\0.975165608\0.975032494\
0.975032494\0.974899313\0.974899313\0.97476595\0.97476595\0.974632414\0.974498892\0.974498892\0.974498892\0.974365187\0.974365187\0.974231447\0.974097694\0.97396386\0.97396386\0.973829968\0.973695705\0.973695705\0.973561302\0.973561302\0.973426812\0.97329233\0.973157792\0.973023165\0.972888502\0.972753686\0.972753686\0.972618599\0.972483401\0.972348192\
0.972212959\0.972212959\0.972077497\0.971942025\0.971806564\0.971806564\0.971671016\0.971535469\0.971399832\0.971263727\0.971127505\0.971127505\0.971127505\0.971127505\0.970990911\0.970854255\0.970854255\0.970716143\0.970716143\0.970577898\0.970439634\0.970439634\0.970301283\0.970301283\0.970162586\0.970023864\0.970023864\0.970023864\0.969885043\0.969746157\
0.969746157\0.969607121\0.969607121\0.969467906\0.969328605\0.969328605\0.969189095\0.969049555\0.96890993\0.968770148\0.968630236\0.968630236\0.968490195\0.968350094\0.968350094\0.968350094\0.968350094\0.968209822\0.968069503\0.968069503\0.967929069\0.967788583\0.967648104\0.967507426\0.967507426\0.967507426\0.967507426\0.967507426\0.967366445\0.967366445\
0.967366445\0.967366445\0.967366445\0.967225072\0.967225072\0.967225072\0.967083346\0.967083346\0.966941484\0.966941484\0.966941484\0.966799431\0.966657323\0.966515221\0.966515221\0.966515221\0.966372675\0.966372675\0.966230095\0.966087497\0.965944844\0.965802025\0.965802025\0.965659107\0.965659107\0.965515789\0.965372381\0.965372381\0.9652289\0.965085397\
0.964941748\0.964941748\0.964941748\0.964941748\0.964797845\0.964797845\0.964797845\0.964797845\0.964797845\0.964652345\0.964652345\0.964652345\0.964652345\0.964652345\0.964506431\0.964360374\0.964214292\0.964214292\0.964067992\0.963921559\0.963921559\0.963774931\0.963628306\0.963481699\0.963481699\0.96333502\0.963188259\0.963041277\0.962894248\0.962746922);
E3=(0.962746922\0.962599438\0.96245194\0.962304394\0.962156773\0.962156773\0.962008591\0.961860378\0.961860378\0.961711958\0.961711958\0.96156345\0.96141486\0.961266259\0.961117297\0.960968325\0.960968325\0.960968325\0.960819143\0.960819143\0.960669959\0.960520666\0.96037136\0.96037136\0.96037136\0.96037136\0.96037136\0.96037136\0.96037136\0.960220987\
0.960070585\0.959919945\0.959919945\0.959919945\0.959919945\0.959919945\0.959768593\0.959617193\0.959617193\0.959617193\0.959617193\0.959617193\0.959617193\0.959617193\0.959617193\0.959464998\0.959464998\0.959311842\0.959311842\0.959158581\0.959005319\0.958851942\0.958698579\0.958698579\0.958698579\0.958698579\0.958545052\0.958391494\0.958391494\0.9582379\
0.9582379\0.9582379\0.9582379\0.9582379\0.958083902\0.958083902\0.958083902\0.957929622\0.957774955\0.957774955\0.957774955\0.957774955\0.957774955\0.957619965\0.957619965\0.957619965\0.957464847\0.957464847\0.957464847\0.957309328\0.957309328\0.957153637\0.957153637\0.957153637\0.956997804\0.956997804\0.956841937\0.956686065\0.956686065\0.956686065\
0.956530084\0.956373803\0.956373803\0.95621749\0.956061106\0.955904677\0.955748076\0.955591451\0.955434795\0.955277962\0.955277962\0.955121\0.955121\0.955121\0.955121\0.955121\0.955121\0.954962774\0.954804098\0.954645324\0.954645324\0.954645324\0.954645324\0.954486493\0.95432756\0.95432756\0.95432756\0.95432756\0.95416855\0.95416855\
0.954009522\0.954009522\0.953850254\0.953850254\0.953690764\0.953531238\0.953371712\0.953371712\0.953212034\0.953212034\0.953212034\0.953052195\0.952892355\0.952732504\0.952732504\0.952572378\0.952412089\0.952412089\0.952412089\0.952412089\0.952251457\0.952090828\0.952090828\0.95192976\0.95192976\0.951768682\0.951768682\0.951607477\0.951446252\0.951284944\
0.951123624\0.951123624\0.950962256\0.950800848\0.950639444\0.950639444\0.950639444\0.950477625\0.950477625\0.950477625\0.950315425\0.950153204\0.949990404\0.949827405\0.949664392\0.949664392\0.949664392\0.949664392\0.94950083\0.94950083\0.94950083\0.949337114\0.949173378\0.949009514\0.949009514\0.94884555\0.94884555\0.94884555\0.948681413\0.948681413\
0.948681413\0.948681413\0.94851707\0.94851707\0.94835196\0.948186731\0.948186731\0.948186731\0.948021228\0.948021228\0.948021228\0.947855307\0.947689386\0.947689386\0.947523355\0.947357331\0.947357331\0.947357331\0.947190969\0.947024487\0.946857953\0.946857953\0.946857953\0.946691255\0.946524514\0.946357651\0.946190793\0.946190793\0.946190793\0.946023761\
0.94585657\0.94585657\0.945689233\0.945521884\0.945354446\0.945186936\0.945186936\0.945186936\0.945019075\0.945019075\0.944850821\0.944850821\0.944682248\0.944682248\0.94451354\0.944344834\0.944344834\0.944344834\0.944344834\0.944344834\0.944344834\0.944175357\0.944005645\0.944005645\0.943835795\0.943835795\0.943665924\0.943495963\0.943495963\0.943325882);
E4=(0.943325882\0.943155626\0.942985392\0.942985392\0.942814967\0.942814967\0.942814967\0.942814967\0.942644357\0.942644357\0.942644357\0.942473429\0.94230248\0.942131355\0.941960133\0.941960133\0.941788768\0.941788768\0.941617217\0.941445654\0.941274104\0.941102559\0.941102559\0.940930621\0.940930621\0.94075829\0.940585684\0.940585684\0.940585684\0.940585684\
0.940585684\0.940585684\0.940412685\0.940239648\0.940066612\0.940066612\0.939893345\0.939893345\0.939893345\0.9397196\0.939545671\0.939545671\0.939545671\0.939545671\0.939545671\0.939545671\0.939545671\0.939545671\0.939370907\0.939370907\0.939370907\0.939196015\0.939196015\0.9390211\0.9390211\0.938845833\0.938845833\0.938845833\0.938845833\0.938670447\
0.938670447\0.938670447\0.938494837\0.938319253\0.938319253\0.938143619\0.93796798\0.93796798\0.93796798\0.937792057\0.937615924\0.937615924\0.937439603\0.937263072\0.937086451\0.937086451\0.937086451\0.936909665\0.936909665\0.936732791\0.936555948\0.936555948\0.936555948\0.936555948\0.936555948\0.936377819\0.936377819\0.936377819\0.936199497\0.936199497\
0.936020453\0.936020453\0.936020453\0.936020453\0.936020453\0.935840668\0.935840668\0.935840668\0.935660462\0.935480231\0.935299883\0.935299883\0.935299883\0.935299883\0.935299883\0.935119065\0.935119065\0.935119065\0.935119065\0.934937369\0.934937369\0.93475556\0.93475556\0.934573728\0.934573728\0.934391818\0.934391818\0.934209879\0.934209879\0.934209879\
0.934209879\0.934027171\0.934027171\0.933844279\0.933661317\0.933661317\0.933478056\0.933478056\0.933294766\0.933111166\0.933111166\0.932926983\0.932742769\0.932558451\0.932374008\0.932374008\0.932374008\0.932374008\0.932188981\0.932003552\0.932003552\0.931818076\0.931632571\0.931447047\0.93126135\0.93126135\0.93126135\0.93126135\0.93126135\0.931075055\
0.930888744\0.930888744\0.930888744\0.930701713\0.930701713\0.930701713\0.930514165\0.930326462\0.930326462\0.930326462\0.930326462\0.930138065\0.929949486\0.929949486\0.929949486\0.929949486\0.929760152\0.929570745\0.929570745\0.929570745\0.929380944\0.929380944\0.929380944\0.929380944\0.929380944\0.929380944\0.92918916\0.928997392\0.928997392\0.928997392\
0.928805459\0.928613126\0.928613126\0.928613126\0.928613126\0.928420021\0.928226946\0.928226946\0.928226946\0.928226946\0.928226946\0.928226946\0.928226946\0.928226946\0.928226946\0.928226946\0.928226946\0.928226946\0.928031824\0.928031824\0.928031824\0.928031824\0.927836316\0.927836316\0.927640467\0.927444576\0.927444576\0.927444576\0.927444576\0.92724812\
0.92724812\0.92724812\0.92724812\0.927051222\0.926853248\0.926655181\0.926655181\0.926456979\0.926258616\0.926258616\0.926060232\0.926060232\0.925861012\0.925661703\0.92546209\0.925262475\0.925062826\0.925062826\0.924862888\0.924662622\0.924462351\0.924261758\0.924261758\0.924261758\0.924060941\0.924060941\0.924060941\0.924060941\0.924060941\0.923859719);
E5=(0.923658372\0.923456846\0.923456846\0.923456846\0.923456846\0.92325444\0.92325444\0.923051693\0.923051693\0.923051693\0.922848018\0.922644281\0.922644281\0.922440225\0.922440225\0.922440225\0.922235278\0.92203034\0.92203034\0.92203034\0.921824808\0.921824808\0.921824808\0.921619012\0.921413184\0.921207316\0.921207316\0.921207316\0.921207316\0.92100099\
0.92100099\0.92100099\0.920794493\0.920587902\0.920381004\0.920174134\0.920174134\0.919966596\0.919758588\0.919550543\0.919550543\0.919342352\0.919134183\0.919134183\0.918925895\0.918925895\0.918717354\0.918508672\0.918299597\0.918090276\0.918090276\0.918090276\0.917880584\0.917880584\0.917670665\0.917670665\0.917670665\0.917670665\0.917460044\0.917460044\
0.917460044\0.917460044\0.917248916\0.917248916\0.917248916\0.917037134\0.916824573\0.916611331\0.916611331\0.916611331\0.916611331\0.916396435\0.916396435\0.916396435\0.916180894\0.91596538\0.91596538\0.91596538\0.91596538\0.915748037\0.915748037\0.915748037\0.915748037\0.915529365\0.915529365\0.915309984\0.915309984\0.915089997\0.915089997\0.914869669\
0.914869669\0.914869669\0.914869669);

J1=(1\0.999877285\0.999754544\0.999754544\0.999754544\0.999631656\0.999631656\0.999508664\0.999385572\0.999385572\0.999262395\0.999262395\0.999139173\0.999015856\0.999015856\0.999015856\0.999015856\0.999015856\0.998892061\0.998768276\0.998768276\0.998644338\0.998644338\0.998644338\0.99852026\0.998396037\0.998396037\0.998396037\0.998396037\0.99827162\
0.99827162\0.998147122\0.998147122\0.998022525\0.997897659\0.997772801\0.997772801\0.997772801\0.997772801\0.997647672\0.997647672\0.997647672\0.997647672\0.997647672\0.997647672\0.997647672\0.997647672\0.997647672\0.997521468\0.997395249\0.997395249\0.997395249\0.997395249\0.997268511\0.997141769\0.997141769\0.997141769\0.997141769\0.997141769\0.997141769\
0.997141769\0.997141769\0.997141769\0.997141769\0.997141769\0.997141769\0.997141769\0.997141769\0.997013857\0.997013857\0.997013857\0.997013857\0.997013857\0.996885746\0.996885746\0.99675757\0.99675757\0.996629046\0.996500415\0.996371794\0.996371794\0.996371794\0.996371794\0.996371794\0.99624301\0.996114192\0.996114192\0.996114192\0.996114192\0.996114192\
0.996114192\0.99598468\0.995855165\0.995725597\0.995595925\0.995595925\0.995595925\0.995595925\0.995595925\0.995595925\0.995595925\0.995465907\0.995465907\0.995335604\0.995205271\0.995205271\0.995074821\0.995074821\0.995074821\0.995074821\0.994944218\0.994944218\0.994944218\0.994944218\0.994944218\0.994944218\0.994944218\0.994944218\0.994944218\0.994944218\
0.994944218\0.994812814\0.994812814\0.994681269\0.994681269\0.994681269\0.994681269\0.994549161\0.994549161\0.994549161\0.99441697\0.99441697\0.99441697\0.994284413\0.994284413\0.994284413\0.994151553\0.994151553\0.994151553\0.994018498\0.994018498\0.994018498\0.993885265\0.993885265\0.993885265\0.993885265\0.993885265\0.993885265\0.993885265\0.993751532\
0.993751532\0.993751532\0.993617235\0.993617235\0.993617235\0.993482859\0.993482859\0.993482859\0.993482859\0.993348291\0.993348291\0.993213656\0.993078862\0.993078862\0.993078862\0.993078862\0.99294386\0.992808446\0.992808446\0.992808446\0.992672703\0.992672703\0.992536714\0.992536714\0.992536714\0.992536714\0.992536714\0.992400502\0.992400502\0.99226416\
0.992127771\0.992127771\0.992127771\0.991991101\0.991991101\0.991854247\0.991854247\0.991854247\0.991854247\0.991716994\0.991716994\0.991716994\0.991716994\0.991579218\0.991441147\0.991441147\0.991441147\0.991302878\0.99116455\0.99116455\0.991026071\0.990887606\0.990887606\0.990748967\0.990748967\0.990748967\0.990748967\0.990748967\0.990748967\0.990610021\
0.990610021\0.990470985\0.990331936\0.990331936\0.990192776\0.990053625\0.990053625\0.9899143\0.9899143\0.9899143\0.9899143\0.9899143\0.9899143\0.9899143\0.989774681\0.989774681\0.989774681\0.989774681\0.989634881\0.989634881\0.98949498\0.98949498\0.98949498\0.98935486\0.989214646\0.989214646\0.989074368\0.989074368\0.989074368\0.989074368);
J2=(0.989074368\0.989074368\0.988933633\0.988933633\0.988933633\0.98879268\0.98879268\0.988651572\0.988510447\0.988510447\0.98836905\0.98836905\0.988227536\0.988227536\0.98808586\0.98808586\0.987943669\0.987801425\0.987801425\0.987658938\0.987658938\0.987658938\0.987658938\0.98751613\0.987373315\0.987373315\0.987230311\0.987230311\0.987230311\0.987230311\
0.987086981\0.987086981\0.987086981\0.986943431\0.986799798\0.986656134\0.986656134\0.986512299\0.986368466\0.986224546\0.986224546\0.986224546\0.986224546\0.986224546\0.986224546\0.986080241\0.986080241\0.985935786\0.985935786\0.985791137\0.985791137\0.985791137\0.985646236\0.985646236\0.985501229\0.985356192\0.985356192\0.985211015\0.985065817\0.985065817\
0.985065817\0.984920047\0.984774054\0.984774054\0.984774054\0.984774054\0.984774054\0.984774054\0.984627712\0.984481332\0.984481332\0.984481332\0.984334645\0.984334645\0.984334645\0.984187717\0.984187717\0.984040565\0.983893359\0.983746108\0.983746108\0.983746108\0.983746108\0.983746108\0.983746108\0.983598124\0.983450134\0.983450134\0.983450134\0.983450134\
0.983301738\0.983301738\0.983153231\0.983153231\0.983004596\0.983004596\0.983004596\0.982855852\0.982707103\0.982707103\0.982558273\0.982558273\0.982558273\0.982558273\0.982409289\0.982409289\0.982409289\0.982260033\0.982260033\0.98211057\0.98211057\0.98211057\0.98211057\0.98211057\0.98211057\0.98211057\0.981960382\0.981960382\0.981960382\0.981960382\
0.981960382\0.9818098\0.9818098\0.9818098\0.9818098\0.981658873\0.981658873\0.981658873\0.981658873\0.981658873\0.981658873\0.98150718\0.981355389\0.981203556\0.981203556\0.981203556\0.981051451\0.981051451\0.980899036\0.980899036\0.980899036\0.980746444\0.980746444\0.980593667\0.980593667\0.980593667\0.980440701\0.980287738\0.980287738\0.980287738\
0.980134532\0.980134532\0.979981207\0.979981207\0.979981207\0.979827611\0.979827611\0.979827611\0.979827611\0.979827611\0.979827611\0.979673456\0.979673456\0.979673456\0.979518884\0.979364324\0.979209767\0.979209767\0.979209767\0.97905505\0.97905505\0.97905505\0.97905505\0.97905505\0.978900127\0.978745173\0.978590153\0.97843509\0.97843509\0.978279837\
0.978124543\0.977969254\0.977813891\0.977813891\0.977658162\0.977502357\0.977502357\0.977346366\0.977346366\0.977190322\0.97703417\0.97703417\0.97703417\0.97703417\0.976877834\0.976721499\0.976721499\0.976564931\0.976564931\0.976564931\0.976564931\0.976564931\0.976407876\0.976407876\0.976250519\0.976250519\0.976250519\0.976092957\0.976092957\0.976092957\
0.976092957\0.975935151\0.975777142\0.975619126\0.975619126\0.975460489\0.975301645\0.975142652\0.974983388\0.974983388\0.974823932\0.974664418\0.974504677\0.974344931\0.974344931\0.974344931\0.974344931\0.974184654\0.974184654\0.974184654\0.97402423\0.97402423\0.97402423\0.97402423\0.973863623\0.973863623\0.973863623\0.973863623\0.973863623\0.973863623);
J3=(0.973702293\0.973702293\0.973702293\0.973702293\0.973702293\0.973540153\0.973540153\0.973540153\0.973377514\0.973377514\0.973214761\0.973214761\0.973214761\0.973214761\0.973214761\0.973214761\0.973051484\0.972888192\0.972888192\0.97272488\0.97272488\0.97272488\0.97272488\0.972561372\0.972397741\0.972233909\0.972070078\0.971906094\0.971741968\0.971741968\
0.971741968\0.971741968\0.97157724\0.97141228\0.971247066\0.971081784\0.971081784\0.971081784\0.9709163\0.970750726\0.970585117\0.970419496\0.97025382\0.970087889\0.969921971\0.969921971\0.969755552\0.969755552\0.969588543\0.969588543\0.969588543\0.969588543\0.969588543\0.969421233\0.969253941\0.96908662\0.96908662\0.96908662\0.968919122\0.968919122\
0.968751569\0.968584016\0.968416359\0.968248694\0.968248694\0.968080887\0.967912949\0.967912949\0.967912949\0.967744501\0.967575948\0.967407365\0.967238764\0.967238764\0.967070012\0.966901269\0.966901269\0.966732415\0.966563563\0.966563563\0.966394448\0.966394448\0.966225159\0.966055746\0.966055746\0.965886288\0.965886288\0.965886288\0.965716729\0.965547174\
0.965547174\0.965547174\0.965377159\0.965377159\0.965377159\0.965377159\0.965377159\0.965377159\0.965377159\0.965377159\0.965206376\0.965206376\0.965035481\0.964864256\0.964692999\0.964521516\0.964350053\0.964350053\0.964350053\0.964350053\0.964177297\0.964004525\0.963831778\0.963831778\0.963831778\0.963658887\0.963485998\0.963313092\0.963313092\0.963140174\
0.963140174\0.962967139\0.962967139\0.962793848\0.962793848\0.962793848\0.962793848\0.96262023\0.96262023\0.962446273\0.962272325\0.962272325\0.962272325\0.962272325\0.962098219\0.962098219\0.962098219\0.961923424\0.96174856\0.961573714\0.961573714\0.961573714\0.961398572\0.961398572\0.961223286\0.961223286\0.961047939\0.961047939\0.961047939\0.961047939\
0.961047939\0.960872253\0.960872253\0.960872253\0.960872253\0.960696398\0.960520282\0.960520282\0.960344039\0.960167734\0.960167734\0.960167734\0.960167734\0.960167734\0.960167734\0.959990343\0.959812589\0.959634839\0.959634839\0.959457073\0.959279325\0.959279325\0.959279325\0.959279325\0.959101175\0.959101175\0.958922981\0.958744797\0.958744797\0.958566468\
0.958388117\0.958209794\0.958209794\0.958030939\0.958030939\0.958030939\0.95785173\0.957672466\0.957672466\0.957492956\0.957313431\0.957313431\0.957313431\0.95713379\0.95713379\0.95713379\0.956954017\0.956774052\0.956774052\0.956774052\0.956774052\0.956593787\0.956413464\0.956413464\0.956413464\0.956413464\0.956413464\0.956232749\0.956052006\0.956052006\
0.956052006\0.95587103\0.95587103\0.95587103\0.95587103\0.95587103\0.955689655\0.955508217\0.955508217\0.955326496\0.955326496\0.955144486\0.955144486\0.954962297\0.954962297\0.954962297\0.954779788\0.954597076\0.954414194\0.954231301\0.954048382\0.954048382\0.954048382\0.953864911\0.953864911\0.953681346\0.953681346\0.953681346\0.953497539\0.953497539);
J4=(0.953313578\0.953313578\0.953313578\0.953129505\0.953129505\0.952945232\0.952760904\0.952576294\0.952576294\0.952391548\0.952206614\0.952206614\0.952206614\0.952206614\0.952206614\0.952021082\0.952021082\0.951835412\0.951835412\0.951835412\0.951835412\0.951835412\0.951649397\0.951649397\0.951462826\0.951462826\0.951462826\0.951275939\0.951089054\0.950902017\
0.950714966\0.95052786\0.95052786\0.95052786\0.95052786\0.950340558\0.950340558\0.950153047\0.949965546\0.949965546\0.949965546\0.949777415\0.949589181\0.949400978\0.949212692\0.949024356\0.948835936\0.948647515\0.948647515\0.948459026\0.948270565\0.948270565\0.948081959\0.948081959\0.947893237\0.947893237\0.947704387\0.947515495\0.947326615\0.947326615\
0.947137713\0.946948593\0.946948593\0.946948593\0.946759248\0.946759248\0.946759248\0.946569773\0.946380199\0.946380199\0.946380199\0.946190252\0.946190252\0.946190252\0.946190252\0.945999896\0.945809399\0.945809399\0.945618786\0.945618786\0.945618786\0.945427935\0.945237011\0.945045991\0.944854867\0.944854867\0.944663384\0.944471836\0.944471836\0.944280201\
0.944280201\0.944088061\0.943895806\0.943703509\0.943511074\0.943511074\0.943318489\0.943125611\0.943125611\0.943125611\0.943125611\0.942932211\0.942738818\0.942545337\0.942351765\0.942351765\0.942158006\0.941963858\0.941769602\0.941769602\0.941575129\0.941575129\0.94138054\0.94138054\0.941185958\0.941185958\0.940991283\0.940991283\0.940796305\0.94060109\
0.940405643\0.940405643\0.94021016\0.94021016\0.94021016\0.940014094\0.940014094\0.939817859\0.939817859\0.939817859\0.9396209\0.9396209\0.9396209\0.9396209\0.9396209\0.939423536\0.939225901\0.939028065\0.939028065\0.939028065\0.938829382\0.938829382\0.938829382\0.938829382\0.938829382\0.938630398\0.938431416\0.938232203\0.938032962\0.938032962\
0.938032962\0.937832993\0.937633039\0.937633039\0.937432695\0.937231967\0.937231967\0.937231967\0.937030775\0.936829613\0.936628459\0.936628459\0.936628459\0.936427016\0.936225589\0.936024068\0.936024068\0.936024068\0.935821855\0.935619655\0.935619655\0.935416967\0.935214117\0.935010744\0.934807377\0.934603907\0.934603907\0.934603907\0.934400281\0.934196621\
0.934196621\0.934196621\0.933992242\0.933787762\0.933582999\0.933582999\0.933582999\0.933378123\0.933173258\0.932968293\0.932763349\0.932558038\0.932352688\0.932147036\0.931941303\0.931735599\0.931529917\0.931324037\0.931324037\0.931118114\0.930912138\0.93070607\0.93070607\0.930499878\0.930499878\0.930499878\0.930293461\0.930086896\0.929880137\0.929880137\
0.929673167\0.929466177\0.929259036\0.929259036\0.929259036\0.929259036\0.929051026\0.929051026\0.929051026\0.928842877\0.928842877\0.928634347\0.928634347\0.928634347\0.928634347\0.928634347\0.928634347\0.928424622\0.928424622\0.928424622\0.928424622\0.928424622\0.928214193\0.928003671\0.928003671\0.927793068\0.92758246\0.927371675\0.927160675\0.927160675);
J5=(0.927160675\0.927160675\0.926948875\0.926736824\0.926524741\0.926524741\0.92631232\0.92631232\0.926099272\0.925885736\0.925885736\0.925885736\0.92567187\0.92567187\0.925457369\0.925242741\0.925242741\0.925242741\0.925027475\0.924812136\0.924812136\0.924596415\0.924380661\0.924380661\0.924380661\0.924380661\0.924164658\0.923948666\0.923732599\0.923732599\
0.923516365\0.923299989\0.923299989\0.923299989\0.923299989\0.923299989\0.923083088\0.923083088\0.923083088\0.923083088\0.922865545\0.922865545\0.922865545\0.922647837\0.922647837\0.92242982\0.92242982\0.92242982\0.92242982\0.92242982\0.922210775\0.9219915\0.9219915\0.921771962\0.921771962\0.921551966\0.921331785\0.921111452\0.921111452\0.920890917\
0.920670431\0.920449623\0.920449623\0.920228342\0.920006364\0.920006364\0.920006364\0.920006364\0.91978216\0.919557671\0.919333012\0.919333012\0.919107533\0.918881934\0.918881934\0.918881934\0.918655995\0.918429711\0.918202833\0.918202833\0.917975265\0.917746666\0.917518009\0.917518009\0.917288907\0.917288907\0.917059209\0.917059209\0.916829168\0.916829168\
0.916598516\0.916367679\0.916136657);

#delimit cr

// assemble pieces of initial data into single vectors

E=E1\E2\E3\E4\E5
J=J1\J2\J3\J4\J5

// Create matrix G (used below) - log{i-1}-log{i}, subset |top_left_row,top_left_column\bottom_right_row,bottom_right_column|
E1=E[|1,.\rows(E)-1,.|]
E2=E[|2,.\rows(E),.|]
G=ln(E1)-ln(E2)

// input initial values of parameters

coeffx=(0.55021\2.83511\1.99822\1.4775\-0.86736\0.70063\0.39241\0.9137)
coeffdx=(0.47666\3.53291\1.43216\0.00704\0.09148\0.97352\0.11888\0.45355)
mxbeta=24.728399807068
mdxbeta=20.5660997877815

// Define new vector containing risk scores
riskvector=J(rows(_inputvars),1,.)
riskvectoropt=J(rows(_inputvars),1,.)
riskvectornor=J(rows(_inputvars),1,.)

// loop over observations and compute risk score for each one

for (i=1; i<=rows(_inputvars); i++) {

	// select observation i of input data
	Xrow=_inputvars[| i,1\i,.|]

	//--------------------- Risk score ---------------------

	// Q - make input vector into column vector
	X=Xrow'

	// R - X*beta
	xbeta=X'*coeffx
	dxbeta=X'*coeffdx

	// $W$5
	expxmx=exp(xbeta-mxbeta)

	// $X$5=EXP(U10-X2)
	expdxmdx=exp(dxbeta-mdxbeta)

	// building vector M - askla{i}==$W$5*(-LN(E2)) and askla{i}=F2*K2*$W$5*G2
	M1=expxmx*(-ln(E[1,1]))

	// F - acs**exp(x-mx)=POWER(E2,$W$5)
	F=E:^expxmx
	// K - ads**exp(dx-mdx)=POWER(J3,$X$5)
	K=J:^expdxmdx

	// M=askla{i}=F2*K2*$W$5*G2
	M2=expxmx*F[|.,.\rows(F)-1,.|]:*K[|.,.\rows(K)-1,.|]:*G
	M=(M1\M2)
	riskscore=colsum(M)
	riskvector[i]=riskscore

	//--------------------- Optimal risk score ---------------------

	// Q - make input vector into column vector
	Xopt=X[1]\X[2]\ln(110)\ln(160)\ln(60)\0\0\0

	// T26 - Xopt*beta
	xoptbeta=Xopt'*coeffx
	dxoptbeta=Xopt'*coeffdx

	// $W$26=EXP(T35-W2)
	expxoptmx=exp(xoptbeta-mxbeta)

	// $X$26=EXP(U35-X2)
	expdxoptmdx=exp(dxoptbeta-mdxbeta)

	// building vector AB - askla{1}_OPT=$W$26*(-LN(E2)) and askla{i}_OPT=Z2*AA2*$W$26*G2
	AB1=expxoptmx*(-ln(E[1,1]))
	// Z - acs**exp(x-mx) OPT=POWER(E2,$W$26)
	Z=E:^expxoptmx
	// AA - ads**exp(dx-mdx) OPT==POWER(J2,$X$26)
	AA=J:^expdxoptmdx
	// G - log{i-1}-log{i}, subset |top_left_row,top_left_column\bottom_right_row,bottom_right_column| - already made
	// AB=askla{i}_OPT=Z2*AA2*$W$26*G2
	AB2=expxoptmx*Z[|.,.\rows(Z)-1,.|]:*AA[|.,.\rows(AA)-1,.|]:*G
	AB=(AB1\AB2)
	riskscoreopt=colsum(AB)
	riskvectoropt[i]=riskscoreopt

	//--------------------- Normal risk score ---------------------

	// Q - make input vector into column vector
	Xnor=X[1]\X[2]\ln(125)\ln(180)\ln(45)\0\0\0

	// T26 - Xnor*beta
	xnorbeta=Xnor'*coeffx
	dxnorbeta=Xnor'*coeffdx

	// $W$39=EXP(T48-W2)
	expxnormx=exp(xnorbeta-mxbeta)

	// $X$39=EXP(U48-X2)
	expdxnormdx=exp(dxnorbeta-mdxbeta)

	// building vector AB - askla{1}_OPT=$W$26*(-LN(E2)) and askla{i}_OPT=Z2*AA2*$W$26*G2
	// building vector AF - askla{1}_NOR=$W$39*(-LN(E2)) and askla{i}_NOR=AD2*AE2*$W$39*G2
	AF1=expxnormx*(-ln(E[1,1]))
	// AD - acs**exp(x-mx) NOR=POWER(E2,$W$39)
	AD=E:^expxnormx
	// AE - ads**exp(dx-mdx) NOR=POWER(J2,$X$39)
	AE=J:^expdxnormdx
	// G - log{i-1}-log{i}, subset |top_left_row,top_left_column\bottom_right_row,bottom_right_column| - already made
	//askla{i}_NOR=AD2*AE2*$W$39*G2
	AF2=expxnormx*AD[|.,.\rows(AD)-1,.|]:*AE[|.,.\rows(AE)-1,.|]:*G
	AF=(AF1\AF2)
	riskscorenor=colsum(AF)
	riskvectornor[i]=riskscorenor
	}

// export risk score vectors into pre-created Stata variables

st_store(., "hriskcvdl30", touse, riskvector)
st_store(., "hriskcvdl30opt", touse, riskvectoropt)
st_store(., "hriskcvdl30nor", touse, riskvectornor)

}

end